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Computation of Influence Coefficients for 
Aircraft Structures with Discontinuities and 
Sweepback 


SAMUEL LEVY* 
National Bureau of Standards 


(1) ABSTRACT 


A general method for computing influence coefficients is out- 
lined which uses Castigliano’s theorem, together with a stress 
analysis of the airplane structure based primarily on equilibrium 
conditions. It is shown, by means of examples, that this method 
is particularly suited to wings having discontinuities, cutouts, and 
sweepback. It is also particularly suited to taking account of end 
fixity at the wing root. 

A method is also presented for use of the influence coefficients 
in computing, by iteration, the normal modes of vibration of an 
airplane as a free body in either symmetric or antisymmetric 
motion. This method avoids the necessity for knowing the loca- 
tion of the elastic axis and gives directly the coupled torsion and 
bending modes. 

Influence coefficients were computed for a sweptback box beam 
(angle of sweepback = 37°) with a large cutout. The influence 
coefficients were nearly twice as large as those for a comparable 
straight beam. 

Computations of relative twist were made for structural units 
consisting of box beams of several lengths in order to indicate the 
optimum ‘size of structural unit for computing influence coeffi- 
cients in torsion. The results indicate that the structural 
unit should have one or two full bays to either side of a large 
cutout. 


(2) INTRODUCTION 


T= PAPER DESCRIBES a portion of a research pro- 
gram on influence coefficients of aircraft structures 
which is being conducted at the National Bureau of 
Standards with the financial assistance and cooperation 
of the Bureau of Aeronautics, Navy Department. The 
purpose of the program is to develop a method of com- 
puting influence coefficients for aircraft structures 
having cutouts and sweepback and to check the method 
experimentally by static and dynamic tests. 





Presented at the Structures Session, Fifteenth Annual Meeting 
LAS., New York, January 28-30, 1947. 
. ° . 
In charge, Dynamic Stress Analysis. 


Only the first part of the program has been com- 
pleted. 
puting influence coefficients of typical aircraft structures 


A general method has been developed for com- 


with discontinuities, such as bulkheads and cutouts, 
with and without sweepback and with a wide range of 
loading conditions. Tests remain to be made to deter- 
mine the range in which the method gives reliable re- 
sults. 

The influence coefficients of an elastic structure are 
defined as the deflections of selected points on the struc- 
ture when a unit force is applied to any one of these 
points. The array or matrix of the influence coeffi- 
cients, therefore, describes the elastic stiffness of the 
structure. 
tribution of mass in the structure, one can compute the 


Knowing this matrix and knowing the dis- 


normal modes of vibration of the structure and their 
frequencies. ! 
used to analyze the behavior of a structure, such as an 
airplane, during landing impact® and during flutter.® 


~4 The normal modes, in turn, may be 


The present report describes a general method of 
computing influence coefficients for stressed-skin struc- 
tures with discontinuities characteristic of aircraft de- 
sign. The method makes use of the stress analysis of 
the aircraft structure based primarily on equilibrium 
considerations,’~'’ together with Castigliano’s energy 
theorem for deriving displacements under load from the 
elastic energy stored in the structure. In those cases 
where equilibrium conditions are not sufficient in them- 
selves to specify the stress distribution, use is made of 
the fact that the actual distribution corresponds to a 
minimum of the strain energy. 

Numerical examples are presented for the case of 
box beams with large cutouts, sweepback, and D-nose 
sections. 
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(3) CASTIGLIANO’S THEOREM 


Castigliano's theorem provides a convenient means 
for computing the deflections of an elastic structure 
under a variety of loading conditions. Castigliano’s 
theorem proceeds from the expression for the elastic 
energy |’ stored in the structure by the application to 
of forces P;, Ps, ..., torques 7), To, ..., 


points 1, 2, ... 
and gives the follow- 


and bending moments J/,, M/s, ... 
ing equations for the displacements ),, 2, along 
P;, Po, ..., the angles of twist a, a2, ... about the 
axes of 7), 72, ..., and the angles of rotation 4;, 6, ... 
about the axes of 14;, Ms, ...: 
\, = OV /OP;, hyo = OV /OP», ce ‘| 
a= OV /OT,, a = OV /OT», ae 
6, = 0V/oM, 6. = OV /0Mz, | 


(1) 


(4) ELastic ENERGY 

The elastic strain energy V is computed from the 
stress distribution in the structure. The elements of an 
airplane structure are, in general, subjected to a plane 
stress distribution. The strain energy stored per unit 
volume with such a stress distribution is 


Vo = (1/2E) (0,2 + ¢,?) — (u/E)o,o0, + (1/2G)r,," (2) 
where 

E = Young’s modulus 

m = Poisson’s ratio 

G = shear modulus, £/2(1 + yz) 

o;, 6, = normal stresses in x, y direction 

Try = Shearing stresses parallel to x, y axes 


By choosing the arbitrary directions x and y in the air- 
plane structural element properly, it is usually possible 


to have o, = 0, giving 
Vo = (0,°/2E) + (t2y*/2G) (2a) 
(5) STRESS DISTRIBUTION 


In an airplane structure of the stressed-skin type, the 
individual sheet areas, spar-caps, and stiffeners have 
little ability to resist lateral forces or torsion. This 
fact makes it possible, in many cases, to determine the 
stress distribution from equilibrium conditions alone. 
Several methods of stress analysis based on this con- 
sideration are presented in references 7 to 12. In 
reference 12, the computed stresses are checked experi- 
mentally for a typical monocoque box in torsion. In 
those cases where the structure is redundant, the condi- 
tion that the strain energy will be a minimum for the 
true stress distribution may have to be used in addition 
to the equilibrium conditions. 

Computations indicate that a stress distribution that 
satisfies equilibrium conditions and differs from the true 
stress distribution by less than about + 10 per cent will 
_ have a strain energy that differs from the true strain 
energy by less than 1 per cent in most cases. The rela- 
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tive accuracy of the strain energy is due to the gradual 
variation in strain energy with stress distribution in 
the neighborhood of the minimum value that corre- 
sponds to the true stress distribution. Castigliano’s 
theorem for the displacements is based on a knowledge 
of the strain energy only, and therefore displacements 
so computed will be considerably more accurate than 
the stress distribution from which the strain energy 
was computed. This fact permits some simplification 
in accounting for redundancies in the structure and in 
determining the effect of the local stress disturbance in 
the neighborhood of cutouts and other discontinuities. 

For purposes of this report we will make the following 
simplifying assumptions in computing the stress dis- 
tribution. It is believed that these assumptions ap- 
proximate actual conditions in many typical stressed- 
skin structures. The degree of approximation remains 
to be explored by tests of actual structures. 

(a) The effective area of the sheet may be lumped 
with that of neighboring reinforcements for purposes of 
calculating axial loads. 

(b) Shear stresses in sheet elements may be com- 
puted from equilibrium considerations. Stresses may 
be assumed uniformly distributed along the edges for 
rectangular sheets and distributed as described in refer- 
ences § and 10 for trapezoidal or warped sheets. 

(c) Changes in axial load in flanges and longitudi- 
nals are the result of the distributed shear in the sheet. 

(d) Bulkheads, spar-caps, and vertical stiffeners 
offer no resistance to lateral forces or to torques. 

(e) Sheet areas offer no resistance to lateral forces. 

(f) Circumferential rings have no flexural rigidity 
in bending perpendicular to their plane and no torsional 
rigidity. 


(6) Loap APPLICATION 


Dynamic loads are applied to the structure because of 
the accelerations of the structure itself or of large 
masses, such as engines, attached to the structure. 

Near the point of application of the dynamic load, 
the structure is primarily subjected to local shear forces. 

At distant parts of the structure, the effect of the 
dynamic load can be resolved into the conventional 
combination of shear force, bending moment, and 
torque. 

(7) Enpb Frxity 

At the root of the airplane wing, free warping of the 
cross section is prevented by the stiffness of the fuselage. 
The effect of this end fixity on the influence coefficients 
for the wing may be taken into account by considering 
a self-equilibrating system of loads applied axially to 
the root end of the wing. The displacements for this 
system of loads are then found by Castigliano’s theorem 
and set equal to zero. This will give equations from 
which the axial loads corresponding to end fixity at the 
root may be computed. 
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COMPUTATION OF 


(S) DEFINITION OF INFLUENCE COEFFICIENT 


The influence coefficient 6,,, between two points m 
and » on an airplane wing is here defined as the dis- 
placement at m due to a unit load at m when the junc- 
tion of the fuselage and wing is considered clamped. 
As a result of Maxwell’s reciprocity theorem, 6,,, = 


0 nm 


(9) STRUCTURAL UNITS 


In order to reduce the work of computing the influ- 
ence coefficients, the airplane will be considered as being 
divided into structural units. The division into struc- 
tural units will be so arranged that, at the cut between 
successive units, the stress distribution is that corre- 
sponding to the conventional beam and shell theory, 
In mak- 
ing this subdivision, it is well to keep in mind the 
argument in reference 9 that local effects may ex- 
tend for great distances in reinforced monocoque struc- 


which neglects the effects of discontinuities. 


tures. 

Fig. 1 illustrates the subdivision of the airplane wing 
into structural units. First, a reference line A—A, Fig. 
1, will be drawn from the junction of the fuselage and 
wing normal to the fuselage. Then planes a, }, c, ..., 
perpendicular to line A—A, will be used to divide the 
wing into structural units I, II, III, .... 

The displacements of points on a given plane a, J, c, 
..., Will then be computed as the cumulative effect of 
the distortion of each of the structural units between 
that plane and the fuselage. Appendix I gives an ex- 
ample to show how the influence coefficients of a com- 
posite structure may be computed from those of its 
structural units. 


(10) Use oF INFLUENCE COEFFICIENTS 


The points m and n, between which influence coef- 
ficients are determined, are the points at which the dis- 
tributed mass of the airplane is assumed to be concen- 
trated. Possible mass concentration points are indi- 
cated by 1M, M2, M3, ... on the sketch of an airplane, 
Fig. 1. Appendix II shows how to choose the location 
and magnitude of these concentrated masses. 

Although the influence coefficients are computed for 

the case of rigid clamping of the wing at the root, they 
may be applied to any other root condition as long as 
the root section remains plane. As an illustration, it 
is shown in Appendix III how these coefficients may be 
used to determine the normal modes of vibration of the 
airplane as a free-free body. These modes may involve 
a combination of torsion and bending and they may be 
symmetrical or antisymmetrical in respect to the y-z 
plane, Fig. 1. 
An obvious advantage of computing modes from 
influence coefficients is that it dispenses with the need of 
determining the elastic axis and gives directly the shape 
of the coupled torsion and bending modes. 
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Fic. 1 


Division of wing into structural units. 


(11) APPLICATIONS 


An application of the foregoing methods to some rela- 
tively simple structural units is given below to show in 
detail how the stress distribution is chosen and how 
Castigliano’s principle is used. In these units, the 
spar-cap area will include the effective sheet area. 


(11.1) Straight Box Beam with Cutout Under One Trans- 
verse Load 

First consider an aluminum-alloy rectangular canti- 
lever box beam with three bays having the bottom sheet 
of the center bay cut out, Fig. 2a. 

The load eonsists of a downward shear force P; on 
the near side of the free end together with a system of 
forces at the root sufficient to hold P,; in equilibrium. 
It is desired to find the deflection at P). 

The loads acting on the various sheet and spar-cap 
elements of the box are determined as follows: Starting 
at the near side of the end, force P;, Fig. 2a, divides into 
two unknown portions a and d in pieces 5 and 1, respec- 
tively, so that 


The remaining shearing forces on the edges of piece 1 
are then determined from the condition that it is in 
translational and rotational equilibrium with the results 
shown in Fig. 2a. At the junction of pieces 1 and 2, no 
external force acts so the forces on these adjoining 
edges must be equal and opposite as shown. The re- 
maining forces on piece 2 are then determined as they 
were for piece 1. This process is continued until the 
edge shearing forces on all the pieces of sheet have been 


assigned. In some cases, such as the junction of pieces 
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Fic. 2. 


thick: bo* cross section, 40 by 10 in.; length of each bay, 40 in. 


Rectangular cantilever box beam with three bay under shear load 


(Bottom of center bay cut out. All sheet, 0.05 in. 


over: ” length, 120 in.; spar-cap area, including the effective sheet 


area, 3 sq.in. per corner; Young’s modulus = 10.5 X 108 lb. per sq.in. 


9, 10, and 14, where it is not certain how the forces are 
distributed, an additional unknown c is introduced as 
shown in Fig. 2a. 

The axial forces in the spar-caps are readily deter- 
mined from the shear forces on the adjacent sheet. 
Consider, for example, the spar-cap at the junction of 
pieces 2 and 3. At the free end of the box, the com- 
pressive force in the spar-cap is zero. The distributed 
shear along the edges of pieces 2 and 3 builds this com- 
pressive force up linearly to (—8) at a distance of one 


bay from the end. Similarly, the compressive force 
in the remaining pieces of spar-cap are determined. 
At the root end, the compressive force on the spar-cap 
at the junction of pieces 12 and 13 will be taken as Ri, 
and similarly for pieces 13 and 14 as Ry, giving 


8b — 8c = R, t (4a) 
12a + 4b + Sc = RJ 


Similarly, we get for the shears at the root end 
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4c = —R; | 
a+b+c= R3/ (4b) 
c= —R;, j 
Combining Eqs. (3), (4a), and (4b) gives 
a= P, —_ (R;/8) o (Rs/4)) 
b = (R,/8) — (R;/4) (5) 
c= — R;/4 \ 


From Eq. (2a), it can be shown that the energy stored 
in a rectangular piece of sheet in equilibrium under 
shearing forces u and v on the edges is 


1.3(uv/ Et) (6) 
where 
u = force on one pair of opposite edges 
v = force on other pair of opposite edges 
t = sheet thickness 
E = Young’s modulus 
0.3 = Poisson’s ratio 


Similarly, it may be shown that, for a piece of spar-cap 
in which the axial force varies linearly from u at one 
end to v at the other end, the strain energy stored is 


(1/6AE)(u? + uv + v?) (7) 
where 
! = length of piece of spar-cap 


A = cross-sectional area 
force on one end 
v = force on other end 
E = Young’s modulus 


=> 
= 


The elastic energy in the box was computed as the 
sum of the energies stored in each of the pieces of the 
box, using Eqs. (6), (7), and (5). This.led to the fol- 
lowing expression for the elastic energy as a function of 
the external load P; and the reactions R,, R;, 


I" = (212.6P,? + 7.78R,? + 31.2R;? — 46.5P,R; + 
61.0P,R; _ 26.1R1Rs)10—® (8) 


The displacements associated with forces R, and R; 
must be zero since the beam is assumed fixed at the 
Hence, we have, from Castigliano’s theorem: 


OV/OR, = (15.6R, — 26.1R; — 46.5P;)10-°=0 


root. 


» (C 
OV/OR; = (—26.1R, + 62.3R; + 61.0P,) 10-* = of 
Solving Eqs. (9) for R, and Rs gives 
R; = 0.900P, (10) 


Inserting Eqs. (10) in Eq. (8) and applying Castigliano’s 
theorem again gives the desired displacement at the 
pont of application of P): 


_oV oO 


iP, 
OP, 


= yp, (135.5P1*)10~* = (271.1P))10-* (11) 
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6. = (271.1) 10~° in. per Ib. (lla) 


11.2) Straight Box Beam with Cutout Under Two Trans- 
verse Loads 
As a second application, we will consider the simul- 
The same rectangular box 
One of the loads will be 
The other load 


taneous action of two loads. 
shown in Fig. 2a will be used. 
taken as P, as in the previous example. 
P, consists of a downward shear force on the far side at 
a distance of one bulkhead from the end, as shown in 
Fig. 2b. It is desired to find the influence coefficients 
for the points of application of P; and P». 

The loads in the different sheet and spar-cap elements 
are determined as they were in the previous case. In 
the notation of Fig. 2b 


a+b=P, 

d = P; | 

8b — 8c + 8d = R, | 
12a + 4b + 8c = R, } (12) 

atb+c=R;| 

ns tin 

—4c = R;} 


Combining the first three and the last of Eqs. (12) 


gives 


a = P, + P. — (R,/8) + (Rs/4) 
—P, + (Ri/8) — (Rs/4) 

c = —R;/4 

a= fi, 


(13) 


The elastic energy in the box was computed from Eqs. 
(6), (7), and (13). This led to 


V = (212.6P,? + 409.6P,? + 7.78R,? + 31.2R;? + 
372.0P,P, — 46.5P,R; + 61.0PiR; — 


104.2P2R, + 190.0P2:R; — 26.1R,R;)10-* (14) 


Again, we set the displacements associated with R; and 
R; at zero, since the root is fixed, or 


OV/OR, = (15.6R; — 26.1R,; — 46.5P; — 
104.2P2) 10-* = 0 


15 
OV/ORs = (—26.1R; + 62.3Rs + 61.0P, + _ 
190.0P2) 10-* = 0 
Solving Eqs. (15) for R; and R; gives 
R; = 0.900P, _ 0.814P% f tian 
From Castigliano’s theorem, the displacements are 
bP + 6y2P 2 = OV /oP, = (0/OP;) (135.5P;? + 
54.5P2? + 74.4P,P.2)10—° = (271.1P, + 74.4P2)10~* (16) 
6o,P1 fe b2P 2 = OV /oP2 = (0/O0P2) (135.5P,? a 
54.5P2? + 74.4P,P.2)10—* = (74.4P; + 109.0P,)10~* 
or the influence coefficients are 
6. = (271.1) 10~* in. per Ib.) 
bio = bo, = (74.4) 10~* in. per Ib. ? (17) 
522 = (109.0) 10~* in. per Ib. 








a2 JOURNAL OF THE AERONAUTICAL SCIENCES OCTOBER, 1947 

















Fic. 3. Sweptback rectangular box beam with three bays under shear load. (Bottom of center bay cut out. All sheet, 0.05 in 
thick; bulkhead, 40 by 10 in.; length of each bay, 50 in.; overall length, 150 in.; spar-cap area, including effective sheet area, 3 sq. 10 


per corner; angle between wing and fuselage, 53.12°; sweepback angle, ¢ = 36.88°; Young’s modulus = 10.5 X 10® Ibs. per sq. in. 


SS Box Beam with Cutout Under Two ously discussed. A plan view of a structural unit of 
the beam is shown in Fig. 3a. The loads will be taken 
as a downward shear force P; on the near side of the 


For the third application, we will consider a box beam 
with sweepback but otherwise similar to the ones previ- end and another downward shear force P2 at the back, 
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COMPUTATION OF 


one bulkhead from the end. It is desired to find the 
influence coefficients for the points of application of 
P, and P». 

The loads in the different sheet and spar-cap elements 
are determined as they were in the previous cases. The 
condition of rotational equilibrium for pieces of sheet 
of parallelogram shape reduces to av = bu, where a and 
u are the length and force, respectively, on one pair of 


edges and } and v on the other pair of edges. In the 
notation of Fig. 3b, 
a+b=P, ° 
d = P, | 
10b — 10c + 10d = R, 
15a + 56 + 10c = R, | (18) 
a+t+b+c=R; 
—-c+d=R, 
—4c¢ = Rs j 


Combining the first three and the last of Eqs. (18) gives 
a = P, + P, — (R,/10) + - 


b = —P, + (R; /10) = (R;/4) (19) 
c= —R;/4 | 
d = P» 


From Eq. (2a), it can be shown that the elastic energy 
stored in a piece of sheet in the shape of a parallelogram 
in equilibrium under shearing forces u and v on the 
edges is 


1.3(uv/Et)(1 + 1.538 tan® @) cos ¢ (20) 


where u, v, E, and ¢t have the same significance as for 
Eq. (6) and @ equals the sweepback angle. Eq. (20) 
reduces to Eq. (6) for the straight wing for which ¢ = 0. 

The actual stress distribution in a piece of sheet in 
the shape of a parallelogram is complicated by shear lag 
effects because the axial component of the force applied 
to the sheet by the bulkheads is gradually diffused to 
the spar-caps. Fortunately, however, since the stress 
distribution used in determining Eq. (20) satisfies equi- 
librium conditions and does not differ too greatly 
from the true stress distribution, the use of Eq. (20) 
should give accurate values of computed influence coef- 
ficients. 

The elastic energy in the box was computed from 
equations (6), (7), and (20), using Eqs. (19) to introduce 
the external loads. This gave 


V = (394.3P,2 + 746.2P2? + 9.33R:2 + 
56.49R;? + 712.6P:P2 — 71.26P:R; + 119.0P,R; — 
154.9P.R; + 347.1P2R; — 38.72R:R;)10-* (21) 


Setting the displacements associated with the root forces 
R, and R; to zero, since the root is considered fixed, 
gives 


OV/OR, = (18.66R; — 38.72R; — 71.26P, — 
154.9P:)10-* = 0 
OV/OR, = (—38.72R; + 112.99R; + 119.0P; +( (22) 


347.1P2)10-* = 0 
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Solving Eqs. (22) for R; and R; gives 


R; = 5.648P; + 6.663P, | 


(29. 
R; = 0.8817P; — 0.7887P» | seid 


Introducing Eqs. (22a) in Eq. (21) and applying Castig- 
liano’s theorem gives for the displacements 


bnP; + bP» = OV /OP, = 

(491.1P; + 143.9P2) 10 .| (23) 

bP; + bP, = OV/OP2 = + 
(143.9P, + 186.8P.) 10-6 


or the influence coefficients 
dn = (491.1)10~* in. per Ib. | 


bo = by (143.9)10~° in. per Ib. 
bx. = (186.8)10~* in. per tb.| 


(24) 


Comparison with Eq. (17) shows that the influence coet- 
ficients were nearly doubled by the sweepback. 


(11.4) Straight Box Beam with Cutout and D-Nose Sec- 
tion Under End Torques 

For the fourth application, we will consider a box 
beam with a D-nose but otherwise similar to the one 
previously discussed. An oblique view of the beam is 
shown in Fig. 4a and an exploded view in Fig. 4b. The 
load will be taken as a torque applied to the two ends. 
At the ends, the stress distribution will be taken equiva- 
lent to that found in a long box without cutouts carry- 
ing the same torque. It is desired to find the twist 
between the two ends of the box. 

The loads acting on the various sheet and spar-cap 
elements of the box are determined as they were in the 
previous cases. However, in applying the conditions 
of equilibrium where D-shaped nose sections are in- 
volved, we must, in addition, make use of the fact that 
a uniformly distributed shear force along an arbitrary 
path ABC, Fig. 5, is statically equivalent to a single 
force R having the following properties: 

(a) Its magnitude is equal to the product of the 
shear force per inch on path ABC times the distance / 
between A and C. 

(b) Its direction is the same as that of a line joining 
A and C. 

(c) Its line of action is at a distance d from the line 
joining A and C, where d is obtained by dividing twice 
the area enclosed between line AC and curve ABC by /. 

As an example, consider the use of this statical equiv- 
alence in applying the condition of rotational equilib- 
rium to the bulkhead shown in Fig. 6a. This bulk- 
head intersects shear webs of the box at AC and DE and 
has a D-nosepiece ABC. The distributed shear on the 
bulkhead along path ABC may be replaced by a stati- 
cally equivalent force v acting in a direction parallel 
with AC and at a distance d; from AC. Both v and d, 
may be determined as described. For rotational equilib- 
rium 


ud, = vd, + wl 
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Fic. 4. D-nose box beam with three bays under torsion load. (Bottom of center bay cut out. All sheet, .05 in. thick; length 


of each bay, 40 in.; overall length, 120 in.; 
10.5 X 108 lbs. per sq.in.) 


In terms of the notation of Fig. 4b, the force distribu- 
tion over the end sections, corresponding to the torque 
T, was taken as 


b + e = 0.000437 
a + f = 0.008627 > (25) 
d—-c= 0.008627 


This force distribution was found in a separate analysis 
of a long box without cutouts. In addition, the axial 
load on the ends of the spar-caps is zero under the given 
loading conditions, and the resultant torque carried by 
the bulkheads is assumed to be zero. This gives 


2d = 3f +g 
7.85a = 40c¢ + 37.85 (26) 
7.85g = 40d + 37.85h 


spar-cap area, including effective sheet area, 3 sq.in. per corner; 


Young’s modulus = 


Combining Eqs. (25) and (26) gives 


1.946a — 0.9465 


c= 
d = 0.00862T + 1.946a — 0.946b 

e = 0.000437 — b ras 
f = 0.00862T — a ‘ 
g 


= —0.00862T + 6.893a — 1.893b 
h = —0.02685T + 12.119a — 2.8930 


From Eq. (2a), it can be shown that the energy stored 
in a D-nose-shaped piece of sheet in equilibrium under 
shearing forces on the edges, Fig. 6b, is given by multi- 
plying the strain energy from Eq. (6) by the ratio L/I 
of the developed length to the distance between spar- 
caps—i.e., 


1.3(uv/Et)(L/l) (28) 
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COMPUTATION OF 


The strain energy stored in the bulkhead, Fig. 6a, is 
computed on the assumption that the shear stress in 
the rectangular portion is «/tl (where ¢ is sheet thickness 
and / bulkhead height) and that in the D-nose portion 
it is v/tl. These stresses were chosen because they are 
correct for a considerable portion of the bulkhead, al- 
though they fail to take into account the bending 
stresses that must be present in the neighborhood of the 
shear web AC. On this basis, the strain energy stored is 


u?\ /d, a ) *) 
‘ - 6: 29 
13(F )(7) + 0.65 (F e td 


The strain energy in the box was computed from 
Eqs. (6), (7), (28), and (29) The forces c to h were 
eliminated by using Eq. (27), and the following expres - 
sion was obtained: 


V = 10-* (0.062447? + 5,844a* + 
745.06? — 3,908ab — 21.883aT + 6.689bT) (30) 


The values of a and 6 corresponding to minimum 
energy were determined from 


01//da = (11,688a — 3,908b — 21.8837)10-* = 01 a, 
dV/db = (—3,908a + 1,490b + 6.6897)10-* = of 
giving 


a = 0.003027, b = 0.00343T (32) 


Substituting Eqs. (32) in Eq. (30) and applying Castig- 
liano’s theorem gives, for the twist a between the two 
ends of the box, 


OV oO 
= =, = => (0.040887?) 10-* = (0.08176T)10-* (33 
a= or ar | ) ( (6T) 33) 
In this analysis the exact values of a and } correspond- 
ing to the correct stress distribution were used. To 
determine the effect of deviations from the exact values, 


let us assume approximate values (error 0.7 and 12.5 


‘per cent, respectively) 


Gappr. = 0.0037, appr. = 0.003T (32a) 


The value of the strain energy obtained from Eq. (30) 
would be 


Vappr. = (0.04098777)10~° 


and the approximate twist, appr, between the two ends 
of the box would be 


Gop, = WV appr. /OT = (0/0T)(0.040987T2)10—° = 


(0.0819747)10-* (33a) 


Comparison of dappr,, Eq. (33a), with a, Eq. (33), shows 
an error of only 0.3 per cent corresponding to the error 
of as much as 12.5 per cent in the stress distribution, 
Eqs. (32) and (32a). 


(12) Optimum Size oF StrucTURAL UNITS 


The structural unit should be chosen small enough to 
reduce computations to a reasonable amount. It 
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Fic. 5. Resultant of forces distributed along a curved path. 
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Fic. 6. Distribution of forces on a D-nose bulkhead and on a 


I>-nose-shaped piece of sheet. 


should be large enough to ‘‘wash out”’ the effect of dis- 
continuities on the load distribution at the ends. The 
optimum size of structural unit represents a compro- 
mise between these opposing requirements. 

For a box beam with cutouts subjected to bending 
loads, the optimum size of structural unit can be deter- 
mined from the well-developed theories for shear lag 
as given in references 12 to 14. 

For a box beam with cutouts subjected to torsion 
loads, the optimum size of structural unit could be de- 
termined from the theory in references 12 and 16. This 
theory is, however, not so well known as the shear lag 
theories. It was decided, therefore, to analyze box 
beams of several sections with a large cutout, in which 
the size is changed by increasing the length of beam to 
either side of the cutout. The effect of size shows itself 
in this case by a change in the relative twist of 
bulkheads to either side of the cutout. In these 
analyses, the spar-cap area includes the effective sheet 
area. 


12.1) Box Beam of Square Section in Torsion 


For the first analysis, consider a series of three square 
box beams of aluminum alloy with three, five, and seven 
bays, having the top sheet of the center bay cut out, 
Fig. 7. The torque loads are applied at the bulkheads 
by shear forces distributed around the box in a manner 
similar to that in an infinitely long box. The resulting 
relative twists are shown in Fig. 7. The relative twists 
@1, @2, a3 Of bulkheads 1 (1, nearest to cutout) and 
Q22, O23 Of bulkheads 2 (2, second bulkhead from cutout) 
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cutout when unit torques are applied at bulkheads 77 for square 
box beams with three, five, and seven bays having top of center 
bay cut out. (All sheet, 0.05 in. thick; box cross section, 10 by 
10 in.; length of each bay, 10 in.; bulkheads, 0.05 in. thick be- 
tween bays; spar-cap area, including effective sheet area, 1 sq.in. 
per corner; Young’s modulus = 10.5 X 106 Ibs. per sq.in.) 
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out when unit torques are applied at bulkheads 77 for 1:2 rec- 
tangular box beams with three, five, and seven bays having top of 
center bay cut out. (All sheet, 0.05 in. thick; box cross section, 
10 by 20 in.; length of each bay, 20 in.; bulkheads, 0.05 in. thick 
between bays; spar-cap area, including effective sheet area, 
1 sqg.in. per corner; Young’s modulus = 10.5 X 108 lbs. per sq.in.) 











rad/Ib in rad/Ib in 
%, ONExIOS a, 0.2ITKIO® 
oi wae 




















rad/ib in rad/Ib in rad/Ib in 
%, 0.109x10° a2, 0.197x10°° X53 0.282x10°° 
A ANT &,5 -203 
ais 121 
Bulkheads 4 3 2 \ { 2 3 4 
' / ' 
at 
rad/tb in rad/Ib in rad/ Ib in rad/Ibin. 
%, O.106x10° x, O.190x10° ,, 0.278x10F a, 0.356x10° 
a, -Ne &,,  -201 Ay, 20! 
X, 120 A, +203 
O, 121 
Fic. 9. Relative twist a;; of bulkheads 77 to either side of cut- 


out when unit torques are applied at bulkheads jj for 1:4 rec- 
tangular box beams with three, five, and seven bays having top 
of center bay cut out. (All sheet, 0.05 in. thick; box cross section, 
10 by 40 in.; length of each bay, 40 in.; bulkheads, 0.05 in. thick 
between bays; spar-cap area, including effective sheet area, 1 
sq.in. per corner; Young’s modulus = 10.5 X 108 Ibs. per sq.in.) 
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out when unit torques are applied at bulkheads jj for a D-nose 
box beam with three and five bays having top of center bay cut 
out. (All sheet, 0.05 in. thick; box cross section, 10 by 40 1.; 
length of each bay, 40 in.; length of D-nose, 20 in.; radius of 
D-nose, 5 in.; bulkheads, 0.05 in. thick between bays; spar-caP 
area, including effective sheet area, 1 sq.in. per corner; Youngs 
modulus = 10.5 X 108 Ibs. per sq.in.) 
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COMPUTATION OF INFLUENCE 


when unit torques are applied to bulkheads 11, 22, 33, 
respectively, are changed less than 3 per cent in going 
from a five-bay box to a seven-bay box. For such a 
box beam and cutout, it would be adequate to limit the 
structural unit size to five bays. 


(12.2) Box Beam of Rectangular (1:2) Section in 
Torsion 


For the second analysis, consider a series of three 
rectangular (1:2) box beams with three, five, and seven 
bays, having the top sheet of the center bay cut out, 
Fig. 8. The relative twists ay, a2, a3, @22, a3 are 
changed less than 1 per cent in going from a five-bay 
box to a seven-bay box. For such a box beam and 
cutout, it would be adequate to limit the structural unit 
size to five bays. 


(12.3) Box Beam of Rectangular (1:4) Section in 
Torsion 


For the third analysis, consider a series of three rec- 
tangular (1:4) box beams, Fig. 9. The relative twists 
1, &2, aNd ay are changed less than 5 per cent in going 
from a three-bay box to a seven-bay box. For such a 
box beam and cutout, a three-bay structural unit would 
be adequate. 


(12.4) Box Beam of D-Nose Section in Torsion 


For the fourth analysis, consider two D-nose multi- 
cellplar box beams with three and five bays, having the 
top cover sheet in the center bay cut out, Fig. 10. The 
torque loads are applied at the bulkheads by shear 
forces distributed around the box in a manner similar 
to that in an infinitely long box. The resulting torsion 
influence coefficients are given in Fig. 10. The relative 
twists ay, a2, and a are changed 4 per, cent or less in 
going from a three-bay box to a five-bay box. For such 
a box beam and cutout, it would be adequate to limit 
the structural unit size to three bays. 


(13) PROGRAM OF FUTURE WoRK TO 
CHECK METHOD 


The following additional work is planned to check, by 
analyses and tests, the adequacy of the method of com- 
puting influence coefficients described in the present 
report: 

(a) Apply the method to the analysis of complete 
airplane wings both of conventional type with cutouts 
and of the sweptback type. 

(b) Make bending and torsion tests on simplified 
monocoque boxes having large cutouts, D-nose sections, 
and sweepback, and compare the resulting bending and 
twisting deformations with theoretical values. 

. (c) Make bending and torsion tests on complete 
airplane wings, and compare the measured bending and 
twisting deformation with the computed influence coef- 
ficients derived as item (a) above. 
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APPENDIX I—INFLUENCE COEFFICIENTS OF STRUCTURE 
CONSISTING OF Two STRUCTURAL UNITS 


Consider a sweptback airplane wing divided into two 
structural units, I, II, as shown in Fig. 11. A set of 
reference coordinates x, y, are chosen with their origin 
fixed in the wing root and with the orientation shown 
in Fig. 11. The planes a, 0, c, at the ends of the struc- 
tural units are taken normal to the x axis. A reference 
line /, parallel to the y axis is chosen in plane } at the 
location of a heavy shear web. The distance of this line 
from the x axis is denoted by z. A similar line is 
chosen in plane c. 

Deformation coefficients are then determined corre- 
sponding to the flexibility of each structural unit. The 
deformation coefficients for structural unit I are ob- 
tained by considering it as clamped at plane a and sub- 
jected at plane 5 to the following loads: 

(a) A load P,, along line /,. The load P, is taken 
positive in the y direction. 

(b) <A bending moment M,, applying forces distri- 
buted as in an infinitely long box. The moment /, is 
taken positive from x to y about the z axis. 

(c) A torque T,, applying forces distributed as in an 
infinitely long box. The torque 7, is taken positive 
from z to y about the x axis. 

The displacements and rotations \,, 9, and a, for 
structural unit I are computed as in Section 3. These 
will be functions of P,, M,, and 7, given by the equa- 
tions 








rh 








Fic. 11. Influence coefficients of a composite structure from 
the influence coefficients of its structural components. 
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Np = MrpP, — AuuM, te MrT, 
6, = O1pP, + OiuM, + Oi7rTy 
a arpP, + amM, + airl, 


The coefficients Arp, Aw, .... Q@1r in these equations 
are the desired deformation coefficients for structural 
unit I where \ denotes displacement of line /, in the y 
direction, 9 denotes rotation of plane } about the z axis 
from x to y, and a denotes rotation of plane 6 about the 
x axis from z to y. 

In the subscripts, I refers to structural unit I, P refers 
to displacement and rotation caused by unit shear load, 
M refers to displacement and rotation caused by unit 
bending moment, and 7 refers to displacement and 
rotation caused by unit torque. 

Deformation coefficients for structural unit II are 
determined by considering plane c subjected to the 
following loads: 

(a) A positive load P, along line /, 

(b) A positive moment M,, applying forces dis- 
tributed as in an infinitely long box. 

(c) A positive torque T,, applying forces distributed 
as in an infinitely long box. 

Plane b is subjected to the following loads: 

(a) A negative load —P, along line J). 

(b) A negative moment — [M, + P(x, — %»)], 
applying forces distributed as in an infinitely long box. 

(c) A negative torque — [T, + P.(s. — 2»)], apply- 
ing forces distributed as in an infinitely long box. 

The resulting deformation coefficients for structural 
unit II will be denoted similarly to those for structural 
unit I. Thus Aqrr is the displacement of line /, in the 
y direction with respect to line /, due to unit torque T 
in structural unit IT. 

The influence coefficients, 5,,, for mass points M/,, 
M2, ..., Me (Fig. 11), of the composite structure are 
computed from the deformation coefficients of the sepa- 
rate structural units as follows: 


614 = 62 = 613 oe O16 = 0 
621 = 609 Te cece 6:6 = 0 
533 = Arp + (23 — 2)Arr + 
(23 — 2) [ar + (23 — %)arr] 
634 = 643 = Arp + (24 — 2)Mir + 
(23 — 2) [ar + (24 — 2)arr] 
535 = O53 = Arp + (X_ — Xp)At + (25 — 2») Arr + 
' 3 — 2)[ar + (x. — X)arm + (25 — 2)arr] 
556 = Oss = Arp - (x, me X»)Atar + (6 ie Zy)Atr + Arp + 
(26 — 2-)Atir + (25 — 2) [arp + (%-. — Xp) arn + 
(26 — 2)arr] + (% — 2.)[anp + (2 — 2-)arnr| + 
(x. — %p) [Arp + (x. — X»)O1m + (% — 2») O17] 
566 = Arp + (x. —_ Xp») Ata + (26 - Zy)Atr + Aue + 
(26 — 2.)Aurr + (26 — 2) [arp + (x, — Xp) orm + 
(z¢ — Z»)arr] + (2 — 2.)[anp + (2 — 2,)anr] + 
(%_ — Xp) [Ore + (%- — Xo) O10 + (26 — 2) O17] 


The influence coefficients, 6,,,, computed as outlined 
are correct not only for the case where the end sections 
remain plane but also when warping is present at either 
plane bd or plane c if plane 0 is taken sufficiently far from 


— 
ta 
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discontinuities in structural units I and II so that the 
warping displacements of the two units match. Twist 
and displacement at plane c due to warping at plane } 
are automatically taken account of by the deformation 
coefficients arrp, Q@t1ry, Q1IT and Alp, Army, Artr- 


APPENDIX II—LOCATION AND MAGNITUDE OF 
CONCENTRATED MASSES 


The location and magnitude of concentrated masses 
approximating the distributed mass of the wing should 
be such that the dynamic loads due to the concentrated 
masses produce nearly the same deflection at distant 
parts of the structure as the distributed mass would have 
done. 

Consider an airplane wing, Fig. 12a, vibrating in a 
normal mode. The maximum deflections y in the nor- 
mal mode throughout the cross-hatched region 5S, Fig. 
12a, can be approximated by the polynomial 


y=Qntoetazt+ dx? +ea2? +... (II, 1) 











Fic. 12. Concentrated masses that produce approximately 
the same deflection at a distant part of the structure as the dis- 
tributed mass would have done. 

The influence coefficient 6 between a point A, distant 
from region S, and a point B, inside of region S, cam 
be approximated by the polynomial 


5 = de + dox + cog + dex? + ee? +... (IL? 


where x, z are the coordinates of B, and dy, hn, @, ds, 
€, ... are coefficients depending on the location of 4. 

If we denote by m the distributed mass per unit area 
in region S and by w the natural frequency in the normal 
mode, the maximum deflection y, at A due to the dy- 


namic loads in region S is given by 
ya = w' fof mys dx dz 


Substituting Eqs. (IJ, 1) and (II, 2) in Eq. (II, 3) and 
expanding gives 


(II, 3 
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va = w [SoS aiaem dx dz + JG S (arb, + ab) X 
mx dx dz + fg JS (ayco + a2c;)mz dx dz + 
SsS (ares + ace, + c1¢2)mz? dx dz + ...] (II, 4) 


Concentrated masses, VM, Ms, M3, ... with coordinates 
Xi, 21} Xa, 22; X3, 23; ... (Fig. 12a) would produce at A 
a maximum deflection, y,4’: 


y. 4 w*(Miyi6) + Mey26o + ws a (II, 5) 


ya’ = w*[Mi(ay + bier + ati + ...)(Q2 + 
box, + C221 + es a + M2(a, + bixe a 
C122 + +») (de + dex, + Coz2 + J+ . 


—_ 
oo 


an 
(II, 6) 


ya’ = w*[ayae(M, + Me + ...) + (ibe + 
d2b,) (Myx, + Moxe + ...) + (ace + 
2C,)(Myz, + Moz, + ...) + (ie, + 
dee, + C1C2)(Myz,? + Moz? + ...) +...) 
(II, 7) 


The deflection y,’, Eq. (II, 7), will be equal to yu, 
Eq. (II, 4), if: 


M,+ Me. + pg a SsJS m dx dz 
Myx; + More + ... = SoS mx dx dz 
Mya, + Moe +... = Sg Sme dx dz 
(II 8) 
Myxyz: + Moxere +... = J { mxz dx dz 
Mya? + Mom? +... = fof mz? dx dz 
Myx? + Mox2? + ... = fg Sf mx? dx dz 
With a finite number of mass points, 14,, Mo, ..., it 
is impossible to satisfy all of Eqs. (II, 8). It was con- 


sidered sufficient to satisfy the first six equations by 
placing the mass points as indicated in Fig. 12b. This 
places four mass points at the ends of the region where 
they will coincide with corresponding mass points from 
the adjacent regions and two mass points in the interior 
of the region. For convenience, it is suggested that the 
region include two structural units and that the two 
mass points in the interior be placed at the junction of 
the two structural units. The masses MM, Mo, ..., Ms 
are chosen to satisfy the six simultaneous linear equa- 
tions: 


M+ Me+...+ Me = SoS m dx dz 
Mix; + Max. + ... + Mexs = Sf mx dx dz 
Mig) + Moz + ... + Mere = Sg [mz dx dz 
Mixy? + Mox2? + ... + Mexe? = Jf mx? dx dz )(II, 9) 
Myxy2, + Mox%e+... + Mexcez6 = 
Ss S mxz dx dz 
Mitr? + Moo? + ... + Meze? = Sg f mz? dx dz 


The first three of Eqs. (II, 9) require the total mass 
and the center of gravity of the six concentrated masses 
to be the same as that of the distributed mass. The 
last three of Eqs. (II, 9) require the moments and 
Product of inertia of the six concentrated masses to be 
the same as that of the distributed mass. 


APPENDIX III—USE oF INFLUENCE COEFFICIENTS IN 
COMPUTING FREE-FREE MopEs oF AIRPLANE WING 
STRUCTURE 


Consider a frame of reference coordinates x, y, 2, fixed 
in space and oriented with respect to the airplane as 
shown in Fig. 1. Let xy, Vn, 2, be the coordinates of 
point mass m. Let the motion of the clamped junction 
of the fuselage and wing be described by A, 0, and a, 
where \ is the upward displacement, @ is the rotation 
from x to y about the z axis, and a is the rotation from 
z to y about the x axis. Consider the fundamental 
mode with a frequency w. More accurate values 
v1", Yo’, ... 9’ of the displacements parallel to the y 
axis may be computed by iteration from assumed values, 


V1, Yo, ... ¥, in terms of the influence coefficients and the 
inertia forces due to vibration by: 
1 } : 
- (y’ —-A- Ox —_ a2;) = Maynbin 
™ n=1 
r 

1 
= On’ — A — Ox2 — am) = > } Mnynden f¢ (III, 1) 
Ww 
° n=1 
. . r 
1 
—=(yr’ — \ — Ox, — az,) = MaYndrn 
pe 

n=1 





In order to determine, A, 8, and a, we compute from 
Eqs. (III, 1) the sums A, B, C, given by: 


r 


1 
) Ma(ym’ — \ — 0Xm — 2m) — = 
w? 
m=1 
) Mn > M,Vnimn = A 
m=1 


n=1 
r 


, 1 
Mumm’ — \ — Xm — Am) a 
@ 
m=1 
> MiXm > M.J.ta, = B 
m=1 n=1 
1 
) Mn2m(m’ — X\ — 0Xm — O8m) — = 
a 


m=1 
> Mn2m > MIJaban = C 


m=1 n=1 


> (III, 2) 





We then make use of the conditions that, if the airplane 
is not being accelerated as a rigid body, 


YY Min’ = 0 


n=1 


> Maxade’ = ‘ (III, 3) 
n=1 


> M,zny,' = 0 


n=1 
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Combining the last two sets of equations, 


r 6 
Fe Mn —~ 2, Mnxn | 
w* vedi 
m=1 “= 
r 
Qa 
- a M n2m A 
w? 
m=1 
r 4 
. 7 : 
canes 6 Mi eu os BE M nXm* 
w- o 
aan’ m= y r (IIT, 4) 
Qa 
- MinXmZm = B 
w 
m=1 
¥ r 
d 7] 
= 9 M n2m aes. M mm xi 
aw” wo 
m=1 m=1 
. 
Qa 9 ] 
a M n2m" = C | 
a } 
m=1 


Since the summations in Eqs. (III, 4) are known from 
the mass distribution and geometry of the airplane, 
d/w?, 0/w*, and a/w? can be evaluated. Knowing these 
values, we can solve Eqs. (III, 1) for the desired ap- 
proximate values y;’/w?, yo’/w?, ... y-’/w?: 


vi : nN Ox az, P 
9 salt » + ; ; + ye + MW nS) n 
w* w* Ww” w* 


n ie } 
, 
Vo NK 6x2  a&e 2 
; 9 - —— = = sie M 1¥ n52n 4 (III, 5) 
w* w~ ie) aw” 
n=1 
: r 
’ | 
Ns XN , OX, | a2, 
= = me ot ers + * + Mi nb; n 
w~ w* ® aw” 
n=1 } 


The ratio y,/(y,’/w*) gives an approximate value of the 
square of the fundamental frequency w*, while the 
ratios (¥;’/w")/(y,’/w?)... give normalized values of 
the deflections in the fundamental mode which are more 
accurate than the originally assumed values. 

The higher modes of vibration may be computed in a 
similar manner provided that the lower modes are re- 
moved from the result after each cycle of iteration, as 
outlined in reference 17. 
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Determination of Turbojet Engine Thrust 


from Tailpipe 


Measurements 


RAYMOND W. PYLE* 


Ryan Aeronautical Company 


SUMMARY 


The engine rating curves furnished by the manufacturer of 
turbojet engines are usually not satisfactory for determination 
of the thrust produced by a particular engine in a particular air- 
plane. These curves do not take into account the effects of 
exhaust nozzle size, tailpipe losses, and, in some cases, inlet 
pressure ratios. For these reasons, it has been found to be more 
satisfactory to determine thrust from instrumentation in the jet 
tailpipe or exhaust nozzle. 

A simplified procedure for the determination of jet-engine 
thrust and airflow has been developed to aid in the reduction of 
data obtained by tailpipe or nozzle instrumentation. Equations 
are derived from which working curves may be plotted to show 
clearly the effect of changes in any variable. These curves are 
especially useful in flight-test-data reduction. 


DEFINITION OF SYMBOLS 


A; = area of exhaust nozzle throat (hot), sq.in. 
= area of tailpipe (hot), sq.in. 

jet thrust, Ibs. 

F, = net or propulsive thrust, lbs. 


i 
on 


P = gas static pressure, lbs. per sq.in. 

P, = ambient (atmospheric) pressure, lbs. per sq.in. 

P; = exhaust nozzle throat (jet) static pressure, Ibs. per 
sq in. 

P, = tailpipe static pressure, lbs. per sq.in. 

Pimp, = tailpipe or exhaust nozzle impact pressure, Ibs. per 
sq.in. 


Piyp = average integrated exhaust nozzle (jet) impact pres- 
sure, lbs. per sq.in. 
T = gac atati ? ° 
= gas static temperature, °R. 


T, = ambient (atmospheric) temperature, °R. 

T; = exhaust nozzle throat (jet) static temperature, °R. 

T, = tailpipe static temperature, °R. 

Timp, = tailpipe or exhaust nozzle impact temperature, °R. 

T na, = tailpipe or exhaust nozzle indicated temperature, °R. 

V_— = velocity of airplane, ft. per sec. 

V; = velocity of jet, ft. per sec. 

% = specific volume of gas in exhaust nozzle throat, cu.ft. 
per Ib. 

V, = velocity in tailpipe, ft. per sec. 

W, = air weight flow, Ibs. per sec. 

W, = gas weight flow, Ibs. per sec. 

Wi = fuel weight flow, lbs. per hour 

R= = perfect gas constant, 53.3 

Y = ratio of specific heats at constant volume and pres- 
sure 


GENERAL THRUST EQUATION 


T2 THRUST DEVELOPED by a turbojet engine has 
been shown to be expressed by: 
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Vi + ASP; — Po) (1) 
where the sum of the first two terms is the gross thrust 
and the last term is known as the ram drag. All terms 
of this equation, except the total weight flow and nozzle 
velocity, may be easily determined from standard in- 
strumentation. 

Weight flow and nozzle velocity can be determined 
by an instrumented rake placed across the tailpipe at 
any suitable position. Two positions are often used: 
at the nozzle or in the tailpipe slightly upstream from 
the nozzle. The nozzle rake has the disadvantage of 
combining, for some flight conditions, subsonic and 
sonic flows which may affect the accuracy of the test 
results. For this reason, the method developed herein 
is based on a tailpipe rake rather than a nozzle rake, 
although equations and graphs may be used for meas- 
urements at either position provided applicable pres- 
sure readings are used. 

A series of total pressures is taken across the tailpipe, 
and each pressure is weighed by the effective percentage 
of area it covers. This average total pressure is then 
used in all computations. 

The velocity of discharge through a nozzle: 


z Ta 2 aoe 
26RTimp( — Jf. -(5 —} | (2)t 
imp. 


If total temperature is measured, 


2ERTime(—™ \i-fh / (Fe) ™ 


V= 
(3) 
If gas temperature is used, 
Time. - TOu/rm (4) 


and Eq. (3) becomes: 


ey 


Or if unshielded temperature pickups are used, as- 
suming the commonly used 60 per cent ram rise in 
unshielded probes: 





VY = 


t Applied Fluid Mechanics, by O’Brien and Hickox; 
revised to agree with the nomenclature used herein. 


Eq. 3,9 
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Tes. = T+0.6(Timp. — 1), T= Tina/[0.4 + 0.6(Pimp./P)~””) (6) 


and 





V = V2gRUr/(y — W) Tina. (Pimp. /P)°~ 97 = 11/104 + 0.6(Pimp./ PI} (7) 
If a different thermocouple factor is known, it may be inserted in Eq. (6) or an interpolation may be made be- 
tween the curves. It should be noted that, by substituting P, or P; for P, V then becomes tailpipe velocity or 


nozzle velocity, respectively. 
Weight flow of gases through nozzle, 


W, = A;P,V;/RT, = W. + Wy 
W, = (A,P,/RT,)V 2gRV y/(y — IV TV (Prme./PS —-1 












































W, = A,P,/V2g/RT,V y/(y — 1)V (Pine/P)°"” — 1 (8) 
Or if Tina, is known. 
A,PV 2¢V ¥/(7 - i Pimp. psidiaadiae Pup. ‘clined 
W, = rn ST ecm — 0.2; —— — 0.4 (9) 
V RT ina. P, P; 
From Eqs. (1), (5), and (8): 
Pry (y-1)/7¥ W.V 
p, = 24,P,(—1) (Pa) a] + aye, - Pa - 7S a0 
7-5 P, g 
For subsonic pressure ratios, P; = Pa, 
F, = 2A P,ly/(y — 1)){(Puve./Pa)'~”” — 1] — (WaV/8) (11) 
At 
Sonic NozzLE CONDITIONS 
Above sonic pressure ratios—jet Mach Number at the nozzle, M; = 1.0: 
M? = 10= [2/(7 - 1) }[(Prup./P)~ ~ ij (12) Fig 
P, = Pimp/((y + are = (13) (P; 
2.1 
Combining Eqs. (10) and (13): es 
2A —1 1 P WV = 
wa aay gaa ae 1|+ A,| 7 — & P.|- (14) is k 
(vy + 1)/2) 2 (iy + 1)/2] g 
Pour /Ps — DI - 2 Prup./Po WV 
F, = 2A,P, 4‘ mp./Pa) [y/(¥ dN 1)/2 \ + A;P, s (Pimp./ = < it “2 WV (15) 
[((y + 1)/2)] Udy + 1)/2]” g 7 
from which, sult 
F, = AyP.([2(Prue/Pa)/{ (Cy + 1)/21”-?} — 1) — (WeV/8) (6 } 2 
= test: 
EQUATIONS FOR NORMAL CONDITIONS 2.1984 stan 
For y = ‘4/3 or 1.333, the average value for the ~ VT. a a 
temperature and pressure range normally encountered, tine Mil ae, Daas nto a typi 
Eqs (4), (5), (7), (8), (9), (11), and (16) become, for Pimp.\"” (Pimp. \""" “ 
velocity and weight flow in the tailpipe: PA, 90.6 P, — ea ror.) — 0.4  (18b) * 
V = 117.175VT,V (Pymp/P,)'4 —1 (17a) iene me while 
2.1984 72 / 1/4 
oy W, = S=— PA y(22") ms (=u (1%) abo 
V = 117.175V Tina. map. / Pp)” ——e ie | 6 ae P, P, : =m 
0.4 + 0.6 (Pimp./Pp)” tet e1 
‘sci ecteabsseapeaig aaa Below pressure ratio for sonic flow P; = re Go 
V = 117.175 Timp.V 1 — [1/(Prur./P,)'* (17c) Pose \¥4 WV of must 
— 7 —_,,— hrust = A,P,8 moo —-1|/-— @ 
W, = (2.1984P,A,/WT,)V (Pip./P,)'*—1 (18a) acceptin ( P, ee 
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Above pressure ratio for sonic flow: 


thrust = A,PA 1.25947 Pip. _ 1) Be (20) 
Py g 

Figs 1 and 2 contain a plot of the several functions of 
(Pimp./P). Fig. 3 is a graph of 117.175 ~/T and 
2.1984/+/7 against measured temperature. Fig. 4 
is a graph of the weight flow per unit area per unit 
pressure for the most commonly used case where Tin. 
is known. 


CONCLUSION FROM TESTS 


This method has been found to give satisfactory re- 
sults if an experimental correction factor, K, is applied 
to the gross thrusts computed by Eq. (19) or (20). 
This correction factor should be determined by ground 
tests of the particular installation in a thrust measuring 
stand. If a weight flow correction is known for the 
particular instrumented ring, it should be used in the 
weight flow equation. It is felt, however, that for 
typical tailpipe instrumentation this factor would be 
very nearly unity. 

The assumption of y = 1.333 for all jet conditions, 
while not strictly true, results in only small errors 
(about one-third of 1 per cent maximum) for the range 
of temperatures and pressures encountered in normal 
let engines. 

Good accuracy in the measurement of pressures 
must be obtained, but temperatures are not critical. 
An error of 10°C. in the tailpipe temperature resulted in 
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an error of net thrust of two-tenths of 1 per cent for a 
typical case investigated. 


SAMPLE CALCULATION 


Alt., 20,000 ft.; I.A.S., 400; OAT = O°F.; T.A.S., 
447 m.p.h.; 656 ft. per sec.; nozzle diameter, 12.4 in.; 
area = 0.8386 sq.ft.; tailpipe temperature = 640°C. 
ind.; airplane static pressure, P, = 13.75 in. Hg = 
972.5 lbs. per sq.ft.; tailpipe static pressure, P, = 23.0 
in. Hg = 1626.8 Ibs. per sq.ft.; P,A; = 815.6; P,A, = 
1789.5; Kr = 0.92. 

Tailpipe total pressure (ref. airplane static) (from 
rake). 
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Pressure Area Pressure Area 
182 0.211 167 0.029 
183 0.114 175 0.057 
182 0.086 177 0.086 
182 0.057 179 0.114 
178 0.029 179 0.211 


160 0.006 


Avg. gc = 179.6 in. H,O 
Ge 13.21 in. Hg 
Pimp. = 13.75 + 13.21 = 26.96 in. Hg 
Pra J P's 1.961 
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Letters to the Editor 


Dear Sir: 

Two recent Letters to the Editor, one by Prof. D. T. Williams 
(JOURNAL OF THE AERONAUTICAL SCIENCES, Vol. 14, No. 4, 
p. 236, April, 1947) and the other by B. Szczeniowski (Jbid., No. 
5, p. 310, May, 1947) have been written with regard to a previous 
Letter of the undersigned which appeared in the JouRNAL, Vol. 
14, No. 1, p. 23, January, 1947. 

In replying to these letters we wish to consider them sepa- 
rately. Regarding the letter of Professor Williams, we should 
like to make the following statements: 

(a) A propeller is not a rocket, and, consequently, formulas 
for rocket propulsion cannot be applied to a propeller and vice 
versa. 

(b) Professor Williams correctly defines propulsive efficiency, 
np, as the ratio of “useful work,” or product of thrust and air 
speed, to “‘energy input’’ to the gases ejected from the propulsor. 
This concept of energy input should alone indicate the fallacy of 
using interchangeably equations for jet propulsion units and 
propellers. The energy input for jet propulsion units, as indi- 
cated in our previous Letter to the Editor, is the ‘‘kinetic energy 
input to the jet from the chemical potential energy and is de- 
pendent on the thermal efficiency conditions whereby the chemi- 
cal potential energy is converted into kinetic energy of the jet,” 
plus the kinetic energy in the fuel due to its being carried in the 
propelled unit at the flight velocity, V. 

The energy input to the propeller equals the change in kinetic 


energy of the air stream as it passes through the propeller. This 
energy input will be: 
Input = (W./2g)(V,? — V?) (1) 
where 
Wa = weight of air flowing through propeller disc per second 


V, = velocity of propeller slipstream, with respect to the 
ground, ft. per sec. 

velocity of the propelled unit with respect to the ground, 
ft. per sec. 


V 


The thrust, F, of the propeller will be: 


F = (Wa/g)(Vp — V) (2) 
The power output will then be 
Output = (W./g)(Vp — V)V (3) 
And the efficiency of the propeller is 
9V 9 
= : (4) 
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From Fig. 1: 


gross thrust/K7A;P, = 1.470 
gross thrust = 0.92 X 1.470 X 815.6 = 1,103 Ibs. 
Pimp./Py = 1.1722 


From Fig. 4 for 
T = 640°C. 
W/P,A, = 0.01111 
W, = 19.9 Ibs. per sec. 


Net thrust (W, assumed negligible) : 


F, = 1,103 — (19.9 X 656/32.2) = 698 lbs. 


Using Eq. (4) for the propulsive efficiency of a rocket, Professor 
Williams obtains efficiencies ranging from 0 to © as V is in- 
creased without limit. However, efficiencies greater than 100 
per cent are without meaning, since both the input and output 
or thrust are negative for forward velocities greater than the 
slipstream velocity. Besides that for V — , n, — 2 and not 
Np ©. 

(c) The various propulsive systems have the following dif- 
ferences in operating characteristics: (1) A rocket uses neither 
the kinetic energy nor the velocity head of the air moving relative 
to the propelled body. Consequently, the value of this energy 
[Zs in the authors’ previous letter, Eq. (5)] in this case is equal 
to zero. (2) A propeller fully utilizes this kinetic energy, and 
it is well known that the design of propellers depends on the utili- 
zation of this energy. ; 

It is easily seen that the propulsive efficiency of a propeller 
may be obtained by means of the general method of attack out- 
lined by the authors in the previous letter. This general method 
is valid for all types of propulsive systems where propulsion is ob- 
tained by reaction. In these analyses it is assumed that the 
coefficient of internal efficiency is equal to 100 per cent. Conse- 
quently, the formulas will give the ideal propulsive efficiency or 
the maximum theoretically obtainable. 

Utilizing the general method of obtaining the propulsive ef- 
ficiency of the propeller gives for the total energy of the air mov- 
ing behind the propeller (with respect to the ground) the expres- 
sion: 


Li = [1/(22)]WeV,? 65) 
Since the energy L; = [1/(2g)]W.V? is introduced from outside 


relative to the body, the total input from the propeller must be 
equal to: 
Li = L; — Ls = [1/(2g)]Wa(V,? — V?) (6) 


This value must be consistent with the value in Eq. (1). 
For the output one obtains: 


Output = ZL, — Ly = Li — Ls — Ly 


where the symbols used denote: 
L; = heat equivalent of the kinetic energy in the air entering 
the propeller, relative to the propeller 
L, = heat equivalent of the kinetic energy in the air leaving 
the propeller which is unavailable for thrust because 
of the forward motion of the propeller 


The unusable kinetic energy, L,, contained in the slipstream ™ 
unit time in this case is equal to: 
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Ly = (Wa/2g)(Vp — V)? (7) 
This equation is well known from the theory of the propeller. 
(See F. W. Durand, Aerodynamic Theory, Vol. IV, Division L; H. 
Glauert, Airplane Propellers.) 
Then output equals: 


L = (Wa/2g)2V(Vp — V) (8) 
and thrust is equal to: 
F = L/V = (Wa/g)(Vp — V) (8a) 
The propulsive efficiency is then 
np = L/L, = 2V/(Vp + V) (8b) 


Hence, using a general approach based on kinetic energy rather 
than momentum, the authors have obtained the ideal Froude- 
Rankine propeller efficiency. 

(3) A ram-jet, turbojet, or any air-breathing jet that brings 
the air essentially to rest at the entrance to the combustion 
chamber does not use the kinetic energy L; but uses the velocity 
head (due to the relative velocity) of the incoming air to improve 
the thermal efficiency of the unit. It cannot use both the 
kinetic energy and the velocity head. Consequently, in the 
total energy balance, this amount must be subtracted as was done 
in Eq. (7) of the authors’ previous letter. 

The difference between various jet systems lies in the internal 
energy processes involving compression and combustion, and 
consequently, in the values of thermal efficiencies. To show this, 
let us compare the expression 


np = [2(v; — V)V]/0;2 (9) 


which is the propulsive efficiency of a jet, as derived in the au- 
thors’ previous letter, to the value obtained by Prof. H. Reissner 
(“Systematic Analysis of Thermal Turbojet Propulsion,’’ Jour- 
NAL OF THE AERONAUTICAL SCIENCES, Vol. 14, No. 4, p. 197, 
April, 1947). 

Assuming in Eq. (2) of Professor Reissner’s paper the value 
Pe = fo and using his expression for Q (p. 200), we obtain from 
his Eqs. (2) and (3) exactly the same formula as Eq. (9) above. 

(4) The so-called ducted fan presents a more complicated 
analysis, since the portion of the air passing through the fan 
utilizes L; as kinetic energy, whereas the portion of the air passing 
through the combustion process utilizes L; as velocity head. 
However, taking the input to the fan as given by Eq. (6) and add- 
ing to it the input to the turbine from Eq. (4) of the previous 
letter gives 





Was (Wat + Wy) Wy : 
pet @ —* (V2 — yy 4 HSA vs iP 

a) 9 ™ a “Ta” ° 
where 


W.y = weight of air passing through the fan, Ibs. per sec. 
Wat = weight of air passing through the turbine, Ibs. per sec. 


Similarly the output will be 


Wet ovv, — V) + 
(2g) ~ . 


Wa W. - 
[Het WO, = (#) v| V (11) 
g g 


If in the tail pipe the jet velocity of the turbine unit, v;,, and 
the jet velocity of the fan unit, V,, have an interchange of momen- 
tum sufficient to give a uniform exit velocity, the jet velocity 
will be V, = vj, = v;, and 





Output = 


— (WasVn) + (War + Wy)oje (12) 
(Wes + War + Wy) m 


The propulsive efficiency will then become 





vj; 


— 21Warlo; — V) + (War + Wyo; - WaVIV 13) 
[Wor(v? — V2) + (War + Wy)0}? + W,V?] 
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If the weight of the fuel, W;, is neglected 
2(Was + War)(vj — V)V 


——— - (14 
[Ways(vj;2 — V2) + Wai;?] 


Np 
The values of propulsive efficiency for the ducted fan will ob- 
viously depend on the ratio of the amount of air passing through 
the ducted fan to the amount of air passing through the turbojet. 
The maximum propulsive efficiency will lie between 50 and 100 
per cent, depending on the value of this ratio. From Eq. (14) 
we obtain Eq. (4) or Eq. (9) for Wa: = 0 or Way = 0, respec- 
tively. 
(d) To calculate the efficiency of a rocket, Professor Williams 
should use the formula for rocket efficiency. In his particular 
case when V = CK, where K = log, (m/m), we have: 


2K 2(V/C) 


eS SS ee (15) 
(K2?+ 1) [(V2/C?) + 1] 


Np 
This efficiency is always less than 1. It is not true that this ef- 
ficiency is independent of the fuel used, as Professor Williams 
claims. It depends at each moment on the mass at that moment. 
This mass is equal to the difference between the initial mass and 
the mass of fuel used from the beginning of the motion till the 
moment considered. The last mass undoubtedly depends on the 
velocity C, which, as is well known, varies considerably for dif- 
ferent fuels. For a prescribed drag variation and performance, a 
thrust of a certain magnitude must be supplied. Since this 
thrust equals the product mC, where m in this case denotes the 
mass of fuel consumed per unit time, selection of a fuel with a 
higher jet velocity, C, will require that less fuel needs to be burned 
per unit time and vice versa. We see that 7, is certaialy not in- 
dependent of the propulsive mechanism as Professor Williams 
has demonstrated in obtaining an erroneous answer by the use of 
the propeller efficiency equation applied to a rocket. 

The discussion presented above should help to clarify some of 
the points raised in Mr. Szczeniowski’s letter, but some further 
discussion is necessary. 

Eq. (1) in his letter is presented as ‘‘the formula for propulsive 
efficiency”’ which ‘‘was established many decades ago.” In 
1865, W. J. N. Rankine and, a few years later, R. E. Froude de- 
veloped the equation for the propulsive efficiency of an actuator 
disc. To date this equation has been correctly applied to the 
theory of the propeller, but it does not apply to rockets and air- 
breathing jets. Eq. (1) differs from their equation, which is 
correct when applied to rockets and jets. Between 1920 and 
1930 the equation for the propulsive efficiency of a rocket in the 
form 


np = 2V0;/(V? + 9;?) 


was derived and presented by Dr. Goddard, in this country, and 
by German and Russian investigators, such as Oberth, Sanger, 
Ziolkovski, Scherschevsky, Noordung, etc. Other references 
may be found on page 4 of Sanger’s book on Rockets. 

This equation, after much discussion in Germany at that time, 
was finally accepted by all of the men working in the field of 
rockets. This is attested to by letters from Drs. Sanger, Prandtl, 
Betz, Seewald, Fuchs, Hopf, and others which were received by 
one of the authors in 1936 (Krzywoblocki). Since all of these 
scientists are still living, as far as is known, and nothing has been 
heard to the contrary, it can be assumed that they still feel the 
above equation to be correct. In addition, Dr. Prandtl and his 
assistants checked this formula by experimentation. (See Kort, 
Raketen mit Strahlappanat, ZFM, 1932.) 

It is difficult to say who first derived the correct formula for 
the propulsive efficiency of air-breathing jets, since it was used 
in classified reports in every country engaged in the last war. 
The authors absolutely do not claim that they are the first to 
present it. The only new item that the authors would like to 
ascribe to themselves is the presentation of a method that 
enables one to find the correct coefficient of efficiency in each 
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TABLE 1 as 
THE CHARACTERISTIC FEATURES OF VARIOUS PROPULSIVE SYSTEMS 
(All results are based on authors’ formulas) 

aa a (A/F = air-fuel ratio) 
Z ———-I tem —___—_—. Rocket Air- Breathing Jet Propeller a, iano Ducted Fan 

eat content input expressed W. - i\fwWet W, W 4 

pg ae oak as a (3 1) 92 | [™s i ‘7 oft & (3) 9;2, for 4 = |) (Vp? v(t) (Vp? — V4) + ee + W; ay 





Kinetic energy of the fuel Lz 


(F2)v: 


Kinetic energy of the incom- 
ing air relative to the pro- W. 
pelled body (or imparted ( *)vs, loss in input 
to the air previously - 2 
rest) 
Velocity head of the incom- 4 Not used Used 
ing air 
Kinetic energy of the jet lh G £0; V)? [a er] v 
| 22 
| 
| 
Input from inside the body | Li + Ls Ty bs: Be 
Output Lili + Le — Lal La + Le Ls 
Thrust (positive) F | Ce (* 7 HL 0; - 
& g 





Fiera ies A 
( 2e yi : 0, for FE “> 0 
~ Ge) 

(vj 


(**)r 


Wy na At 
(3) V2 20 for —_ @ 


Was ‘ Wat\,,. 
Ct) - Ger 


Used for turbine only 


| G)v 2 from 
out 
side 

Not used 


(vj — V)2? + 


nl(F2)v, — Vv)? [4 + Wy 
<8 
for —> = Cet Drm 


if vj = 


Ii Lig + Lit + Leet 
Lig + Lit + Le — Ley — La — Lu 


vj and =. > o 


Vp =2 


Ls I1 — Lh 
: dh edna 
(#8), _ (22) ne pe ty 
g g 4 
; W. : 
(er ~ ( =! (ti — V). for Fa 


| 
aa mg ecg de np, neg- | 2V0; 2[(Wa + Wy)vj — WaV]V 2(vj — VV 2V 
ecting the differences in 2 2) | VW, 4 > ™ 72 ‘ = Zas 7 d 
pressure of the undisturbed I Wa + Wee a al , a ed ae or ae oe 
flow at the entrance and for ——_ 
the exit | | 
Limits of np (never negative | < < o< <1 7 
for a positive thrust) slide adi 4 < A <0.5 oa pt. 0 "Pp 1 OS s1 
Maxi , » V¥ , , 
Maximum np | 1 for V = 0; ™ oe ~~ 1 = 4 -§ Ad se 0 1 for Vp = JV 0.5 for War = 0; 1 for Wat = 0 
| P his tia : 
Thrust at maximum np | (Wy/g)V (Wa/2)V, for A/F — @ 0 eee 
(Wy/g)V, for A/F — 0 ° 


propulsive system. This coefficient is different for a propeller, 
rocket, air-breathing jet, ducted fan, etc. 

Mr. Szczeniowski’s reference to the Otto cycle efficiency as 
approaching a limit of 100 per cent is quite true for the air stand- 
ard analysis, which is entirely hypothetical. Similarly, reference 
is made to Carnot’s cycle efficiency as giving a limiting value of 
100 per cent for 7; infinitely large or T, = 0. With the actual 
sink temperature of the atmosphere and the existing limits of 
attainable combustion temperature, the Carnot cycle efficiency 
actually states that any cyclic process which will operate must 
have an actual thermal efficiency well below 100 per cent. 

Since thermal efficiency equations give asymptotic values of 
100 per cent, it does not follow that propulsive efficiency equa- 
tions should give a similar result. Mr. Szczeniowski, without 
proof to the contrary, questions the correctness of ascribing the 
input kinetic energy of the air to the thermal efficiency of an air- 
breathing jet. A clarification of this point may be made by con- 
sidering that the air entering the jet must be accelerated from 
rest, with respect to the ground, to the velocity V of the pro- 
pelled unit, since it was assumed to be brought to a negligible 
velocity with respect to the combustion chamber. The energy 
involved, (W./2g) V2, is imparted to the air by the propelled 
unit. Obviously, this energy, Ls, is a loss as far as the propulsive 
efficiency is concerned. 

A further argument in support of ascribing this energy to the 
thermal efficiency may be had from consideration of the ram-jet. 
If the normal axial air intake were to be replaced by a nearly 
radial air intake, the effect on the thermal efficiency is obvious, 
since the compression ratio would drop nearly to 1 and the jet 
velocity would be very low. 

Regarding the ‘‘best power,”’ it 
breathing jet, which, from the propulsive efficiency standpoint, 
has maximum power when V aad maximum propulsive 
efficiency when V = 1;/2, will also have the optimum value of 
the useful power when V 

The heat content of a pound of hydrogen as against that for a 
pound of gasoline is a matter that should be discussed in con- 
nection with thermal efficiency. 


seems obvious that an air- 
= ) » 


= v;/2 


= 9;/2 


= 0; 


In Mr. Szcezeniowski’s analysis of the energy release by com- 
bustion, he states that sufficient energy is released so that the 
velocity of the fuel relative to the earth is changed from — V to 
w. This is perfectly true, but the energy released by combustion 
is not the total energy output of the system. To see this, take 
the energy of a quantity of fuel before combustion as 

Prota = [W,/(2g)] V? + Q 
where Q is the useful energy released by combustion. The total 
kinetic energy of the products of combustion will be 

KEouwur = [W;/(2g)] W? 

The change in kinetic energy will be due to a portion of the heat 
release duriag combustion and will have the value 

AKE = [W;/(2g)|(w? — V?) 

released during combustion will also ap- 

This will be 


A portion of the energy 
pear as thrust power. 


Prnrust - (W; g)ujV 
The energy Q supplied by the fuel is then 


Q = AKE + Pibrust = [Wy/ 


with w = 7; — V.. The total energy in the system is then 
: Wie. We... - 
Protal = V2?+0 = — (V? + 2;?) 
(29 (20) 
(2g (<2 
and the losses are 
: Wy ~ W, 
Prota = Penn ie a, time V)? = uw 
(2g) (2g) 


Actually, if the fuel were scooped up from rest with respect t 
the ground, the energy (w;/2g) V? would have to be imparted to 
it by the propelled body. This would necessarily reduce the 
energy available to accelerate the products of combustion from 
the velocity V to the velocity w 


(Continued on page 582) 
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Aerodynamic Performance of Delta Wings 
at Supersonic Speeds 


A. E. PUCKETT* ano H. J. STEWART? 
California Institute of Technology 


ABSTRACT 


In recently published papers the authors have discussed the 
problem of drag of delta wings with straight trailing edges‘ and 
also the lift of this type of wing. The present paper consists of 
an extension of these results to delta wings with swept-back or 
swept-forward trailing edges. Drag at zero lift for two slope 
surfaces and lift is computed for these airfoils, and a discussion is 
given of the drag due to lift. Centers of pressure are also com- 
puted. In addition, some results for drag of wings with three 
slope surfaces are given. -Discussion is made of optimum design 
from the standpoint of leading and trailing edge angles, and loca- 
tion of maximum thickness. It is shown that considerable im- 
provements over two-dimensional airfoil performance can be ob- 
tained. Results are presented in the form of conventional engi- 


neering coefficients. 


1. INTRODUCTION 


. HAS BEEN SHOWN by various authors that the super- 
sonic flow past thin or slender bodies can be com- 
puted to a first approximation by assuming that the 
disturbance velocities produced by the body are small 
relative to the free-stream velocity. The equations of 
motion become linear if second and higher powers of 
the disturbance velocities are neglected, and it is thus 
possible to use superpositions of simple solutions to 
obtain more complicated solutions. The linearized 
theory of the two-dimensional wing at supersonic 
speeds was given by Ackeret! and the theory for axially 
symmetric bodies was given by von Karman and 
Moore.’ A third special case, that of conical flow, was 
first investigated by Busemann.* 

One of the authors developed solutions for flow about 
thin wings at zero lift using a distribution of sources in 
a plane, with strength proportional to the local slope 
of the airfoil in the flow direction. The elementary 
solutions used in this case were conical solutions. The 
method was applied particularly to the computation of 
the drag of wings of triangular plan form, or “delta’’ 
wings, with a straight trailing edge. Similar source 
solutions were also found by Jones.° The problem of 
the lift of such a wing was investigated by the other 
author,* using the conical properties of the flow to reduce 
the solution to a conformal transformation problem. 
The limiting case of the lift problem for extremely low 
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supersonic Mach Numbers was first investigated by 
Jones.’ 

It is the purpose of the present paper to generalize 
the previous results to the case of a delta wing 
with swept-back or swept-forward trailing edges. Re- 
sults are given for the delta wing with a double-wedge 
profile. Partial results are also given for the delta 
wing with a three-slope profile. The problems of zero- 
lift drag, lift, center of pressure, and lift-drag ratio are 
treated for the entire range for which the flow is conical. 


2. DRAG oF DELTA WINGS 


In reference 4 the fundamental linearized equations 
for the computation of drag of thin airfoils at supersonic 
speeds were presented and final results were given for the 
delta wing with a straight trailing edge. The general 
equations will not be repeated here, but certain of the 
elementary solutions may be used directly in the com- 
putation of the drag of a delta wing with swept-back 
trailing edge and double-wedge profile, as seen in Fig. 1. 

The initial step in the solution is the computation of 
the pressure distribution over a single delta of constant 
slope, A, as in Fig. 2. This boundary condition is satis- 
fied by a distribution of sources of constant strength 
within the delta. Since solutions to the linearized equa- 
tions may be added, several such deltas may be super- 
imposed to obtain the desired distribution of slopes 
over an airfoil plan form. The pressure distribution 
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‘ VM? —1 = B, n= s/B 
| 7 tel The pressure coefficient was then shown in reference 
4 to be a function only of 6, m, and ¢, given by the real 
t= CONSTANT part of the function 
+r In? = 2 
oe = Rl Ps = cosh 1 Pe (1) 
re x Bev n* — 1 l — & 
This form of the function is convenient for m > 1 and 
t < 1, i.e., for pressures within the delta, when the Mach 
wave is ahead of the delta. When the leading edge is 
NJ ahead of the Mach wave (n < 1), the real part of Eq. 
Fic. 2. (1) is written more conveniently as 
over the wing of Fig. | is composed of superimposed 4d * n2— fp? 
ee, aos ; = a ee ae i Ri| — — sin- @/ - 2) 
pressure distributions associated with each of the deltas p BrV1 mae 5) ear (2) 


having vertexes at A, B, and C. 

It was shown in reference 4 that the pressure over the For the region n < ¢ < 1, that is, within the delta, but 
single delta of Fig. 2 is constant along rays through the ahead of the Mach wave, the only real term in Eq. (2) 
If Cartesian coordinates (x, y) are intro- is the first term in the brackets. This indicates that 


vertex, O. 
duced, then the pressure will be constant along lines = _in this region the pressure is constant, and has, as a 
constant, where ¢ = sy/x and s is the tangent of the matter of fact, the value corresponding to the two- 


dimensional flow normal to the leading edge. 
In the case that m > 1 and >t > 1—that is, in the 
region behind the Mach wave and ahead of the leading 


sweepback angle, ¢. 
The pressure coefficient may be defined as 


C= (b — Pr)/C/s)aU" edge of the delta—the real part of Eq. (1) is 
where /;, pi, and U are pressure, density, and speed, a . 
respectively, in the undisturbed stream. Let d be the C=— eo y ioe. (3) 
7 BrV n? — | t? — | 


slope of the surface of the delta in Fig. 2, and further let 


The basic pressure distributions are now available in Eqs. (1), (2), and (3). The contribution to the drag of the 


pressure on any surface element dA, will be the pressure times the local slope, dD/qg = C,\dA, where gq is the dy- 
namic pressure. It must be remembered that the \ appearing in this equation is the actual local slope of the air- 
foil, while the \ appearing in the expression for C, is associated with the source strength of the particular delta 


The integral of the basic pressure distributions of Eqs. (1) and (2) over areas of 


whose effect is being considered. 
Since the elementary 


the type shown in Fig. 3 will be required in the computation of the drag of the complete wing. 
area of Fig. 3 is 


dA = A[(1 — r)?/(1 — tr)? |dt 


where A = c*/2s, we must compute 


4X(1 — 7)? : ] n* — ft? 4A ; 
f C,dA = — 2 af ———. cosh? _ - dt = G(n, r) (4) 
Ai BrvV n? — 1 o (1 — #r)? l—F TB 


Here all functions of » and 7 are collected in G, which may be evaluated as 


Stivd ‘1-7 logn 4 r cosh—! n 4 2 ' , V1 — rn? (5 

7(",7r) = = —— tan ——————_> 
l+r] Vn? —1 Vn? — 1 V1 — rn? n—m+Vn? — 1 

This form is most convenient, of course, form > 1 and rn < 1. 

venient form for the real part of this expression is: 


2V rn? ma ] 


‘i —? log n r cosh! n 
G(n, r) = 


] 
_- - ~ aaa + — lo 1 + ——€ 5 ———— 
l+r]Vn2-1 Vn? —1 Vrn? — 1 | n(l —r) + Vn? —1 —V Pn? — 1 


and in the case that » < 1, the real part becomes 


l-—r r l 7 : 7 
ty =_ - ——SES —1 “ om is —1 7) 
G(n, r) = “gr pth PS apa cos~! nm + a 3 + sin im | ( 





In the case that m > 1 but rn > 1, a more cot- 
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Since the source distribution in the deltas with vertexes at B and C, Fig. 1, may influence the pressures in regions 
ahead of their leading edges, the integral of the pressure distribution of Eq. (3) over an area of the type shown in 
eference Fig. 4 is required. If the pressure coefficient of Eq. (3) is written in terms of the m defined by the leading edge of 

the real the area in Fig. 4, we have 


4r l — 7 2 rn l 24,2 a l 
' / Ch = peti f ——. cosh-!4/——— 
(1) As rBV rn? — eer how 


where B = rc*/2s. This may be written as 
> 1 and 





he Mach Sas C,dA => 4NBF(n, r)/rB (8) 
edge is where 
t of Eq. ; , re 
: l-—r log rn l rn? — 1 + V (r?n? — 1)(m? — 1) {1 
ee eck SA. a a A th «lant (9) 
2 l+r Vn — I Vn? — 1 n(1 — r) \ 
) (2) 
eIta, but The drag of the complete double-wedge delta wing of 
Eq. (2) Fig. 1 may now be composed of proper combinations Tet  T-A-wING DRAG 
tes that of the functions in Eqs. (4) and (8). The total drag will TT | (MAX THICKNESS POSITION 
as, aS a consist of contributions in regions (1) and (2) due to v 0-0 
he two- the source distributions in deltas A, B, and C. Wemay 
let Cp,4 = drag in region (1) due to source distribution 
s, in the in A, with a similar notation for the other terms. If 
leading the slopes in regions (1) and (2) and the source- 


strengths in A, B, and C are expressed in terms of a, r, 
and r, which describe the trailing edge sweepback, loca- 
(3) tion of maximum thickness, and thickness ratio (see 
Fig. 1), we finally have: 


277 (1 —a 
ig of the Crs = heed G(n, r) 
pad! af (1 —a) 
the ai. me or] A 2 lee a) — G(n, r)] 
ar delta BrL(l — r)(r — a) 
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aaa ae — 7)(r — a)? a(m, ) 
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_ 2r) atl —a@) st a 
Coo = ae i? a 5 to ae P (m2) 


at 27° a(1 — a) a 
Corre  v | = 2 lF(™2) ° 


The sum of these six terms gives the drag coefficient 
for both surfaces of a symmetrical double-wedge airfoil 
with the plan form shown in Fig. 1, at zero lift, based 
on the actual plan area of the wing. 

It is also convenient to describe the location of the 
maximum thickness point in terms of fraction of root 
chord, using 


Coz = 





6 = (1 — r)/(1 — a) 


It is then clear that the drag coefficients for this fam- 
ily of delta wings at any Mach Number may be ex- 
pressed most generally as 


(CpV M? — 1)/r? = f(n, a, b) 


In order to reduce the number of interesting inde- 
pendent variables, it is worth while first to plot the 
drag coefficient for given plan form, i.e, given values of 
n and a, as a function of }, the location of maximum 
thickness. This is done in Figs. 5 and 6 for a = 0 and 
0.5. It will be noticed in both cases that a minimum is 
reached for values of } in the vicinity of 0.15 to 0.20, 
independently of the value of m. Further computa- 
tions for a = 0.25 and 0.75 have shown this to be sub- 
stantially true in those cases also. Thus an optimum 
configuration will be obtained with the maximum thick- 
ness at 15 to 20 per cent of the root chord for any double- 
wedge delta wing. The discontinuities in the curves of 
Figs. 5 and 6 occur when rn = 1, i.e., when the maxi- 
mum thickness line is parallel to the Mach wave. 
When an is greater than 1, this, of course, will not 
occur. , 

The conditions under which this location of the peak 
point may or may not be important are perhaps shown 
more clearly by plotting Cp@/r? against 1/n for fixed 
values of b. Since 


k = 1/n = (VM? — 1)/s 


this abscissa represents a distorted Mach Number 
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ZA- WING DRAG 
SWEPT TRAILING EDGE 





+ 


TWO DIMENSIONAL DOUBLE WEDGE 


Fic. 8. 


scale for a given value of leading edge sweepback. 
These plots are shown in Figs. 7 and 8 for a = 0 
and 0.5. 

It will be remembered that the drag of a two-dimen- 
sional supersonic airfoil, of double-wedge profile with 
maximum thickness at the mid-chord point, is, accord- 
ing to the linearized theory, 


(CpV M? — 1)/r? = 4 


From Figs. 7 and 8 one observes that at low values of 
k, ie., large sweepback or low Mach Number, the wave 
drag is reduced so much that the value of b does not be- 
come critically important. However, the value of b does 
become critically important as k approaches 1, i.e., the 
Mach wave approaches the leading edge, and it is deci- 
sive in the determination of the largest k at which the 
delta wing is superior to the two-dimensional double 
wedge. 

The family of curves for a = 0.5 (Fig. 8) shows three 
discontinuities in slope: (1) at an = 1, when the Mach 


wave is parallel to the trailing edge; (2) atrm = 1, when 
the Mach wave parallels the maximum thickness line; 
and (3) atm = 1. The last two constitute the impor- 
tant peaks in the curves. 

The effect of sweeping the trailing edge may now be 
studied more easily by plotting curves of (Cp X 
V AM? — 1)/r? vs. 1/n for various values of a, using a 
constant value of b, say b = 0.2. 
plotted in Fig. 9. It is observed there that the effect 
of sweeping the trailing edge is somewhat similar to the 
effect of moving the peak point forward, in that the 
drag at low values of k is reduced and the range in which 
the delta wing is superior to the two-dimensional 


These curves are 


double-wedge wing is extended. 

Another interesting effect is observed in the approach 
of the curves to the point k = 0, corresponding to Mach 
Number 1. For the case a = 0, it was mentioned in 
reference 4 that the drag coefficient for extremely large 
values of m is approximately 


CyV M? — 1/7? = 2 log n/xr?n 
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O DIMENSIONAL DOUBLE WEDGE | __ 


Fic. 9. 


The derivative with respect to k is, therefore, logarithmically infinite at k = 0, and the drag coefficient, for a 


given sweepback, approaches infinity as M approaches 1. 


However, for a greater than zero it can be shown that 


the behavior of the drag coefficient for extremely large values of m is given by 





CoV M?-—-1 21 ] (1 — a) 
—_——_—_ = -- log r(1 — a) — 
7? anLa(l — r)(r — a) 


The derivative with respect to k (= 1/n) of this func- 
tion remains finite and, therefore, the drag approaches 
a finite value as M approaches 1. 

Finally, it will be of interest to plot actual drag coef- 
ficients as functions of Mach Number for several air- 
foils of given geometry. These curves can be easily 
computed from the curves of Figs. 7, 8, or 9. In Fig. 10 
are plotted drag coefficients for wings with 30° and 15° 
semivertex angle (60° and 75° sweepback), each for 
values of a = 0 and 0.5. 

The curves for zero trailing-edge sweepback (a = 0) 
exhibit the two characteristic peaks at the points where 
the Mach waves are parallel to the maximum thickness 
line and leading edge. The weak approach to infinity 
near M = 1isalsoseen. For a greater than zero, the 
peaks are increased in height but the drag at low Mach 
Numbers is considerably decreased. In addition, the 


1 
———— log a(1 


r(r — a) 


1 
— r) * a =i log (r —_ | — a)’ 





singularity at the point where the Mach waves are 
parallel to the trailing edge is also apparent, with an 
appreciable reduction in drag for Mach Numbers below 
this point. 

The foregoing computations have all concerned delta 
wings with a double-wedge profile. It might be 
guessed, from analogy with the two-dimensional 
theory, that this configuration will probably offer the 
least drag for a given thickness ratio. However, from 
a practical point of view, the optimum profile will be 
that shape giving least drag for a given structural 
strength, presumably as measured by the section modu- 
lus of the wing root section. If the structure is con- 
centrated in a single spar at the maximum thickness 
point, then the double-wedge design should be nearly 
optimum. However, another extreme might be con- 
sidered—a solid airfoil, for which the least drag with 












































it, for a 
wn that 


| — a)? 


ves are 
vith an 
s below 


d delta 
sht be 
nsional 
Ter the 
r, from 
will be 
uctural 
modu- 
is con- 
ickness 
nearly 
e col- 
g with 


AERODYNAMIC PERFORMANCE OF 


Fic. 


given section modulus will probably result from a more 
distributed thickness. To obtain an idea of the trend 
of this shape, a single family of profiles is considered, 
with surface composed of three slopes, the central one 
being zero. (See Fig. 11.) The drag coefficient of this 
airfoil will again be composed of six parts, which may be 
written in terms of the functions F and G, and 7 and a, 
similar to the manner in which the drag of the delta 
wing with swept trailing edge was expressed. In Fig. 
12 are plotted curves of (Cp\/ M2? — 1)/r? vs. a, the 
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position of the second maximum thickness line, for sev- 
eral values of r, and for m = 2.5. It is observed that as 
a decreases from r the drag rises, but not rapidly. 

The section modulus for this modified double wedge 
may be written as 


S/C8r? = (1 + 3r — 3a)/24 
We can therefore form the ratio 
W = CpB/(S/C*) 


which is independent of 7, and attempt to find a mini- 
mum value for W, corresponding to the least drag for 
a given S. Curves of W vs. a may be drawn for each 
of the values of r used in Fig. 12; a sample curve for r 
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n=25 
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= 0.9 (and n = 2.5) is given in Fig. 13. It will be ob- 
served that a minimum occurs in the vicinity of a = 0.6, 
i.e., with the second peak at the 40 per cent chord point. 
By plotting similar curves for other values of r, and n = 
2.5, it will be discovered that the minimum value of IV” 
occurs in the vicinity of r = 0.9, i.e., forward maximum 
thickness point at 10 per cent chord. It is observed 
from Fig. 13 that a saving of approximately 17 per cent 
in drag for given section modulus is realized from the 
case of the simple double-wedge profile (@ = r). 

Thus, for a wing of solid construction, the region of 
maximum thickness should be distributed between ap- 
proximately the 10 per cent chord point and the 40 per 
cent chord point, permitting a slightly thinner section 
and reducing the drag. Consideration of a more general 
family of thickness distributions would undoubtedly 
lead to still better profiles but the trend is clearly indi- 
cated by this first step. For practical design purposes, 
it would appear reasonable to round off the corners of 
the modified double-wedge section, and probably to use 
a rounded leading edge, since the local conditions on the 
leading edge (with m > 1) are essentially characteristic 


of subsonic flow. 


3. Lirt oF DELTA WINGS 


The lift of a delta wing with a straight trailing edge 
was calculated by Stewart® for the range 0 < k < | 
where k = +/M? — 1 tan w and by Puckett! for the 
range 1 < k. These results will be expanded in the 
present section to cover the more general plan form with 
swept-back or swept-forward trailing edges over the 
entire Mach Number range for which the flow is conical. 
For this range the pressure distributions computed for 
the straight trailing edge in the earlier papers is still 
applicable. This restriction to conical flows implies 
that the present solution to the lift problem holds for all 
Mach Numbers high enough so that the fact that there 
is a trailing edge does not influence the lift distribution 
over the wing. This limit is (see Fig. 1) 


1< V M? — 1 tan o| (10) 
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For the lower portion of this Mach Number range or, 
more exactly, for the range k < 1, the wing is inside the 
Mach cone through the leading edge point, and the lift 
distribution [see Eq. (45) of reference 6] is given by 


Ap 4 tan? wo 
== E(k’). F <0 SEU ( l 1) 
aq (kV tan? w — tan? w 
where 

Ap = pressure below wing minus pressure above 
wing 

a = angle of attack of wing 

q = dynamic pressure 

k’ = V/1 — FF 


E(k’) = complete elliptic integral of the second kind 


Half of this lift is, of course, the result of an increase in 
pressure below the wing (for a > 0) and the other half 
the result of pressure decrease above the wing. For the 
higher Mach Number range, 1 < k, the Mach cone from 
the leading edge point cuts across the wing and the cor- 
responding lift distribution [see Eq. (35) of reference 4 
or Eq. (2)] is 


Ap + tan w 
aq Vie 1 
2 (M? — 1)— — tan? w 
Rii 1 — — sin - x wa 


T tan? w) — tan? w 
It should be noted that only the real part of the expres- 
sion in brackets in Eq. 12 is significant. 

If the trailing edge of the wing is more nearly normal 
to the mean flow than the leading edge, i.e., if io > wy, 
pressure distributions of both types will exist. If 
io < wo (this can occur only with an exaggerated 
sweep forward of the trailing edge), only pressure dis- 
tributions of the type given by Eq. (12) exist in the 
conical flow regime. The latter case, however, is not 
believed to be of any technical importance. 

If the lift coefficient C;, is based in the usual manner 
on the dynamic pressure and the wing area, the lift 
curve slope, dC,/da, is given by the mean value of 
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Ap/aq [from Eq. (11) or (12)] over the wing area. 
Thus, for k < 1, 


dC, 4 tan? wo ‘wo dA 


da ~ -E(k’)(A/2) + Wee wo — tan? ° 








where the wing area A is given by 
A = ?* {sin o sin w/sin (¢ — w»)] 


and the triangular area element (see Fig. 3) is 








dA = '/,[/? sin? o/sin? (¢ — w) |dw (13) 
From this 
dCy _ 4 tan? wo sin (o — wo )sin o 
penne . x< 10 is 20 25 30 35 40 45 80 $5 60 
da E(k’ ) sin Wo MACH NUMBER 
wo dw Fic. 15. Lift of several A-wings. 
c> oan aa aie 
f sin? (¢ — w) V tan? wy — tan? w ) 


The definite integral may be conveniently evaluated by It may be noted that for a straight trailing edge, ¢ = 
means of the substitution tan w = tan w sin 6. It is 90°,@ =0,and H(a) = 1. Eq. (15) then reduces to the 
then found that result given in Eq. (47), reference 6. It may also be 
: _. noted that the parameter a is purely geometrical. The 

dC,/da = [2m tan wo/ E(k’) \/1(a) (15) minimum Mach Number for which the present theory 

where applies is given [see Eq. (10)] as 1 < k/|al, ie., Eq. 
(15) is valid for the range |a|] < k < 1. For the lower 
Mach Number range, 0 < k < |a| the flow is not of a 
al Ds ies 1 —a = (—a) | conical type and the solution for this case is beyond the 

scope of the present paper. 


-——-—— cos 
l1+a_ (1 —a?)3” 


a = tan w)/tan o 


H(a) 





For the higher Mach Number range, 1 < k, the slope of the lift curve, which is given by replacing Eq. (11) by 
Eq. (12) in Eq. (14), is 


dC, 4sin (o - _ — wo) sin ¢ wo 2. (Mm? — 1-1 — tan? w dw : 
— = . Rl 1 — -sin— - — (16) 
da +/&* — loos wo * cs tan? w) — tan? w sin? (¢ — w) 


This definite integral may be most conveniently evaluated by making the transformation tan w = ¢ tan w. The 
final result, valid for the range 1 < &, is 





da (1 +a)VM?—1 








+a 
V1 — a2/k? V1 — 1/k? 


The final values of dC,/da as given by Eqs. (15) and (17) are seen to be functions of the three parameters k, a, 
and M, and a complete representation would require a family of charts. It may be noted, however, 
that 1/,./M? — 1 (dC,/da) is a function of only the two parameters k and a and may thus be presented on only 
one chart. These results are presented in Fig. 14. The ordinate was taken as '/,»/.M? — 1(dC,/da) because the 
linearized theory of a two-dimensional wing at supersonic speeds gives unity for this quantity. It is of interest 
to note that the delta wing with a sweptback trailing edge has a higher value of dC,/da than the two-dimensional 
wing for a wide range of Mach Numbers. For two particular plan forms, namely: a = 0 and a = 0.5, these re- 
sults are also plotted for various values of wo in the conventional engineering form, dC,/da vs. M, in Fig. 15. As 
it appears that the most important range technically is the range |a! < k < 1, it is useful to rewrite Eq. (15) in 
the following form: 


dC, 8 Ee (— a/k) cos~! A (17) 
é 


(\/ M? — 1/4)(dC,/da) = [xk/2E(k’)|H(a) = J(k)H(a) 


The right-hand side of this equation is seen to be a product of two functions, one containing only a and the other 
containing only k. These functions are tabulated in Table 1 and their product is plotted in Fig. 14. 
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In order to calculate the stability and trim characteristics for a delta wing, it is necessary to know the position 
of the center of pressure of the lift forces. For the range of Mach Numbers considered here, the flow is conical 
and the center of pressure may be readily conputed since the pressure is constant on any radial line through the 
origin, and the center of pressure of any infinitesimal area element as given by Eq. (13) is two-thirds of the length 
of the element from the leading edge. In the present case it is only necessary to compute the center of pressure 
for the range 1 < k, since the center of pressure for the range la| < k < 1is independent of k and equal to the 
value for k = 1 for the range 1 < k since both Eqs. (11) and (12) give the same lift distribution for k = 1. From 
this the pitching moment M of the lift forces about a spanwise axis through the leading edge point of the delta wing 


is given by [see Eqs. (13) and (16)]: 


M 8/8 sin* o tan ay we 2. (M2 i 1)- tan? w| cos wdw 
—= == Ri l — -sn— ae ate aeell iS em (18) 
aq 3V rR? — 1 o | T tan? w — tan?w |sin?(o — w) 








It is most convenient to present this result as a dimensionless center-of-pressure position 


= M/L = (M/aq)/(dC,/da)lA 
i.e., is the number of root chord lengths that the center of pressure is back of the nose. If Eq. (18) is evaluated 
it is thus seen that 
2 if =: eos ( /k) 
ce -—a 
. 1 k? a(3 — a”) cos —!(1/k) ak 
e as , ; ; 
3 (1 — a*)(1 — a?®/k?)*”? (1 — a*)(1 — 1/k*)*/* ~=—k* — a? 
[cos~'! (— a/k)/VW] — a?/k?] + @ [cos (1/k)/V1 — 1/k?] 


(19 















































- This solution is valid for the range | < k. For the 
range la| < k <1, the corresponding result is obtained 
30 by letting k — 1 in Eq. (19). It is 
33 cos~! (—a) 
2" | - a Se a(4 — a*) + (2 + a’)-—>G . — 
| | eal scien a ae 
4 ; a cos~! (—a) 
2 a1 —- &)ia-+ —— 
'0lo— 3 V1 — a 
Ole oe 1 l l 
| It may be noted from both Eqs. (19) and (19a) that for 
0 :, a = 0, & = 7/3. The center of pressure positions for 
5 2.0 2.5 3-0 y I 
° = - ; ; , various values of a and & are presented in Fig. 16, and 
ke Vu? tan w, for |a| < k < 1, #is given in Table 1. 


A third item, in addition to the lift and the center of 
pressure, which is of interest in any wing design is the 
drag due to lift. In the subsonic region this drag 1s 
conveniently called the induced drag and arises from 
the finite span effects. At supersonic speeds this drag 


Fic. 16. Center of pressure of A-wing. 





TABLE 1 
Lift and Center of Pressure of Delta Wing for |a| < k < 1 








Py tea a " iy 4 — is more complex, with one part (sometimes) arising from 
oe ~ f00e ~O4 ee oe 
0.05 0.07816 0.7 114994 2 6168 finite span effects and one part arising from the energy 
0.10 0.15459 0.6 1.3399 2.0219 sropagated in the nose and trailing-edge waves. These 
0.15 0.2284 0.5 12386 1.6713 oe : denies 
0.20 0.2991 0.4 1.1653 1.4418 different parts cannot be clearly separated and will be 
0.25 0.3662 0.3 1.1097 1.2807 spoken of, here, as the drag due to lift. 
0.30 0.4298 0.2 1.0656 1.1621 re asta 
0.35 0.4897 0.1 1.0297 1.0714 As the lift is produced by purely normal pressures, tt 
0.40 0.5460 0 1.0000 1.0000 appears that a thin wing at an angle of attack a from 
0.45 0.5989 —0.1 0.9748 0.9426 PP ‘ ‘ sia 8 . 8 . 
0.50 0.6485 =§).9 0.9531 0.8955 zero lift which has a lift coefficient C, will also have at 
0.55 0.6949 —0.3 0.9340 0.8563 associated drag which is given by 
0.60 0.7384 —0.4 0.9179 0.8232 8 8} y 
0.65 0.7792 -—0.5 0.9029 0.7950 ~ ms o/ Dy) 
0.70 0.8172 -0.6 0. 8887 0.7706 Cp = aC, = C,*/(dC;/da) s 
0.75 0.8528 —0.7 0.8777 0.7495 tad ; —. 
0.80 0.8861 —0.8 0.8670 0.7309 This result is known to apply to the two-dimensiona 
pg pipet supersonic wing theory and also applies to the delta 
0.95 0.9743 wing for the case k 2 1. The proper solution for the 
ubdesed 4.0080 range |a| < k < 1, is not, however, quite so simple. 
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AERODYNAMIC 


The difficulty here depends on the fact that the pressure 
distributions that produce the lift [see Eq. (11)] have 
a singularity at the leading edge. Exactly the same 
result is obtained in the subsonic theory of the flat plate 
airfoil. If the singularity is overlooked, one predicts 
(since the pressure acts normal to the airfoil surface) 
that the lift is accompanied by a drag as given by Eq. 
(20); however, a more careful treatment of this prob- 
lem shows that the drag due to lift is zero. The effect of 
the singularity at the nose is to produce a finite suction 
force that just cancels the drag contribution from the 
rest of the wing. It must thus be expected that the 
singularities at the leading edges of the delta wing for 
the range |a| < k < 1 will produce a suction that will 
reduce the drag below the value calcu!ated by Eq (20) 
No entirely satisfactory theory for the magnitude of 
this effect has been developed; however, an approxi- 
mate theory based on an analogy with the flat plate two- 
dimensional airfoil at zero Mach Number will be pre- 
sented. 

For the flat plate airfoil of chord ¢ the suction force 
per unit span is 


T = xpcU*a? 
Near the leading edge (on top) 
p = —ap V?V/c/x 


where p is the pressure change from the ambient pres- 
sure and x is the distance from the leading edge. As T 
varies with a? while the p varies linearly with a, we 
may take as a general formula for the suction at = 0 


T = (x/pU’)|pr/x]? (21) 


the quantity in brackets measuring the strength of the 
pressure singularity 

The pressure near the leading edge of the delta wing is 
[see Eq. (11)] 


ap U? tan®/? wo 
v 
2 sec wo (k’) 





‘R/x (22) 


where x is again the distance from the leading edge and 
R is the distance from the leading-edge point. If Eq. 
(21) is applied to the local flow normal to the leading 
edge and the Prandtl-Glauert correction is made, the 
suction force is then 

m/1 — M,? 


aU [pvV/x]}? = 


T= 
ra’pU?+/1 — M,? tan wo 
2E*(k’) 


where the U,, and M,, are the velocity and Mach Num- 
ber of the flow normal to the leading edge. The corre- 
sponding decrease in drag* is ' 


6D = — 2 f-® T sin wdR 





(23) 


oasteeenitteneeenesiies 


* This decrease in drag is in agreement with the result pre- 
sented by von Karm4n in the Wright Brothers Lecture for 1946. 
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Fic. 17. Drag due to lift of A-wing. 
or 
ra’pU? sin w tan wV/1 — M,?_, 
6D = — ———— ao —* R? (24) 
2E*(k’) 


If this resulting suction is expressed in dimensionless 
form and subtracted from Eq. (20), it is seen that 


Cp = aC, {1 — [k’/2(1 — a)H(a)E(k’)]} (25) 


This corrected formula for the drag due to lift should be 
used in place of Eq. (20) for the range ja! < k < 1. 
This result is shown in Fig. 17. For a = 0, k = 0.7, 
the effect of the suction is to decrease the drag due to 
lift by 26 per cent. 


4. DISCUSSION 


From the standpoint of drag reduction, it is evident 
that there are two applications in which the delta wing 
may offer advantages. First, and simplest, occurs 
when a wing with fairly small sweepback of both leading 
and trailing edge is used. In this case the strong drag 
peak produced by a two-dimensional wing at Mach 
Number | disappears and is replaced by a weaker peak 
at a higher Mach Number, corresponding to’coincidence 
of the Mach wave with the leading edge. In this ap- 
plication it is presumed that the design Mach Number 
is well above the Mach Number corresponding to this 
peak. Inspection of Fig. 7 or 8 indicates that the 
maximum thickness point should be located near the 
mid-chord point, to reduce the height of the peak, and 
the value of Cp/r? at design Mach Number will be nearly 
equal to the drag of a two-dimensional wing. In this 
application, corresponding to comparatively small 
sweepbacks of 30° to 45° for design Mach Numbers of 
1.5 to 2.0, the primary advantage obtained from the 
delta wing is the reduction in the transonic drag peak. 
It should be remembered, also, that the dC,/da for this 
design, corresponding to k = 1.5 or more, is equal to, 
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TABLE 2 
Comparison of Delta and Rectangular Wings at M = 1.50 














Case 1 2 3 
Delta, Delta 
wo = 20 wo = 20° 
b= 0.2 b = 0.2, 
Rectangular k = 0.4 k = 0.4 
Wing Plan Form A.R, = 1 a=Q0 a = 0.25 
VM? — 1dCr 0.554 0.544 0.591 
4 da 
du 1.98 1.94 2.11 
da 
Relative Area, S 1.00 1.018 0.938 
Relative root chord, 1 1.00 1.69 1.26 
Root thickness ratio, r 0.10 0.059 0.080 
a — 
Cov _ as 4/,* 1.54 1.240 
- 
Cp, (wave drag) 0.0119 0.0048 0.0070 
Co, (skin friction) 0.0060 0.0060 0.0060 
Coy 0.0179 0.0108 0.0130 
SCop, (relative drag) 0.0179 0.0110 0.0122 
(Cr/Cp) max. 5.25 8.6 9.3 





* This value is correct for A.R. > ©. The proper value for 


A.R. = 1 has not been computed. 


or greater than, the two-dimensional value (depending 
on the trailing edge sweepback) and, therefore, the lower 
aspect ratio permits a much larger root chord and lower 
per cent thickness. The reduction in thickness will 
generally be sufficient to reduce the drag of the delta 
wing below that of a rectangular wing at the design 
Mach Number. 

The more spectacular demonstration of the properties 
of the delta wing occurs if the sweepback is sufficient 
that the Mach wave is always well ahead of the leading 
edge, i.e., k = 0.4 to 0.7 at the design Mach Number. 
In this case, inspection of Fig. 14 shows that the value 
of dC,/da, if the trailing edge is swept-back, may be 
nearly equal to the two-dimensional value. On the 
other hand, inspection of Fig. 8 indicates that the value 
of Cp/r? has dropped well below the two-dimensional 
value, while the lower thickness ratio will cause the 
wave drag coefficient to become extremely small. 

These characteristics can be demonstrated by a com- 
parative example of several wings designed for a Mach 
Number of 1.5. Three cases will be considered: (1) a 
rectangular wing of aspect ratio 1, of double-wedge pro- 
file, tapering from a root thickness of 0.10 to 0 at the 
tip; (2) a delta wing with b = 0.2,k = 0.4, anda = 0; 
(3) a delta wing with b = 0.2, k = 0.4, anda = 0.25. 
The areas were adjusted so that the product of dC,/da 
and area is the same in the three cases. This implies 
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that the stabilizing effectiveness of the three wings 
would be the same, or that the lift would be the same 
for a given angle of attack. The thickness ratios of the 
wings were adjusted so that the absolute thickness is 
the same in the three cases. The results of this com- 
putation are summarized in Table 2. 

It will be observed that the wave drag coefficient of 
the delta wings is of the same order as the skin friction, 
and that the final drag coefficients, including an allow- 
ance for skin friction, are roughly two-thirds of the 
minimum drag of the rectangular wing. The maximum 
lift-drag ratio of the delta wings is also seen to be 70 
per cent greater than that of the rectangular wing. 
Since wing weights would be proportional to wing area 
in the first approximation, it is seen that the wing 
weights of these three designs are approximately equal. 

The example given above gives a qualitative indica- 
tion of the relative merits of these wings. It is clear 
that the optimum wing plan form in a specific case can 
only be determined by a complete engineering analysis 
of the design, in which such factors as the type of struc- 
ture, the spanwise center of pressure location, and the 
importance of structural weight must all be given proper 
consideration. 


5. CONCLUSIONS 


Methods have been developed which permit the 
computation of the lift and drag of delta wings for a 
large family of plan forms and thickness distributions. 
It appears in a simple example that the delta wing 
properly designed can be used to reduce the minimum 
drag to roughly two-thirds that of a low-aspect-ratio 
rectangular wing, without impairing its lifting efficiency. 
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The Associated Matrices of Bending and 
Coupled Bending-Torsion Vibrations 


WALTER P. TARGOFF* 
The Glenn L. Martin Company 


SUMMARY 


The associated matrix of the bending vibration of a beam is 
developed. This is combined with the torsional matrix to obtain 
the associated matrix of the coupled vibration. By continued 
premultiplications of the associated matrices of the successive 
beam sections a final matrix equation is obtained relating the 
boundary conditions on both ends of the vibrating beam. The 
imbalance of any chosen end condition is plotted against the 
frequency constant A in fashion similar to Holzer’s method. A 
numerical example is given. Some advantages of the method 
are pointed out. The extension of the analysis to the flutter 


problem is noted. 


INTRODUCTION 


XN THE PRESENT TIME there are in general use three 
methods for the computation of bending and 
coupled bending-torsion natural frequencies of beams. 
The first is the matrix iteration method,! the second (for 
coupled modes) is the classical method of using 
the uncoupled modes as generalized coordinates in 
Lagrange’s equation,” and the third is the method of 
Myklestad.* 

A major objection to the first two methods is that 
they require that all the lower modes be determined be- 
fore a particular desired solution can be obtained. They 
also require completely separate treatment for the sym- 
metric and antisymmetric cases of the airplane wing. 
Myklestad avoids these troubles but introduces com- 
plexities into his solution which are not easily followed 
by the average engineer. 

An analysis of the problem leads to the conclusion 
that trouble in the coupled mode calculation arises 
from the handling of the bending vibration. This 
paper develops a method in which bending frequencies 
can be determined to permit facile coupling with the 
torsional oscillations. 


SYMBOLS 


W = weight of section of beam, Ibs. 

1 = length of section of beam between stations, in. 

I = bending moment of inertia of area of section of beam, 
in.4 

E = bending modulus of elasticity of beam, Ibs. per 
sq.in. 

S = static moment of section of beam, Ib. in. 

y = relative bending deflection of beam, in. 

n = frequency of vibration, cycles per sec. 
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g = acceleration of gravity, in. per sec.? 

AK = frequency factor = 49r2n?/g = 0.1023n? 

|’ = shear force, lbs. 

M = bending moment, lb. in. 

@ = bending slope angle, rad. 

T = torsional moment, Ib. in. 

@ = relative torsional angular deflection, rad. 

J = polar moment of inertia of beam section, Ib. in.? 
C = torsional rigidity of beam section, Ib. in. per rad. 


A, B, D, F = matrices, defined in text 


BENDING CALCULATIONS 


The beam is divided in the usual manner into sections 
with discrete masses lumped at each station, having 
weightless elastic shafting in between (Fig. 1). 

Consider a mass Wy, detached from the beam, vi- 
brating through an amplitude yy at frequency 
(Fig. 2). 

If the shear just to the right of the mass is Vy’, then 
the shear just to the left is Vy’ + KWyyy = Vy. In 
fact, we may write 


Vy = Vy! + KWwyw 


My = My’ 
ox = oy’ 
Vv = Vn’ 


In matrix notation these equations may be written 
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FIGURE | 
EQUIVALENT VIBRATING BENDING SYSTEM 
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FIGURE 2 
CHANGE IN DYNAMIC CHARACTERISTICS ACROSS WEIGHTLESS IWAFT 
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Vy I 0 0 KWy 
My 0 1 0 0 
N = x 
oy 0 0 1 0 
Yn 0 0 0 l 
or [By] = [Ay] X [By’]. 
3). 


Next, consider a beam length /y detached from the beam (Fig. ‘ 
Here, 


TOBER, 1947 
Vy" 
M,’ 
: (1) 
On 
yy’ 
Vy 
x My (2 
4) 
on 
Yn 





Vv+1 = Vy 

M'y +1 = My + Vyly 

oy +1 = Oy + (Myly/Ely) + (Vyly?/2ETy) 

vy +1) = In + Only + (Myly?/2EIy) + (Vrly'/6ETy) 

Again, in matrix form: 

Va +3 1 0 0 0 
M'y +1 a ly 1 0 0 
o'y + 1 ra by? 2EIy ly/EIn l 0 
yn +1 ly? /6EIy ly*/2EIy by l 


or [By + 1] = [Dy] x [By]. 


Combining the above two matrix equations by premultiplying Eq. (1) by Eq. (2) we get, for the complete beam 


section (Fig. 4): 














V'w + 3 I l 
M'y +1 ly 
ly? 
@nt+1 = oR ly 
, ly’ 
r+ 
Snail : | 6ELy 
or 
[B’y + 1] = [Dy] x [By] ad [Dy] x 
[B’y + 1] as [Fyw’] x [By] 
Therefore 


[By’ + 2] = [Fv +1] X [B’y + 1] 


and 
[B sete _ bead [F] xX [Bright onal 


and 


[Brett _ = [Arete — x [F] x [Bright eat] 


atid [Frotai] x [B’ right etl 


In words, by successive premultiplication of the dy- 
namic elastic matrices we arrive at a matrix equation 


relating the boundary conditions at 
beam. 


If our choice of the frequency constant K was 
correct, these end conditions are satisfied. By plotting 
the imbalance of any desired end condition against K 














0 0 KWy 7 Poe 
1 0 KWyly My’ 
ly KWyly? , (3) 
Ely ; m1" 1* 
ly? KWyly® yy’ 
Z ) i 2s st 
2Ely y T 6EIy § : 








the various frequencies can easily be determined as in 
Holzer’s method‘ for torsional systems. 

Obviously, the natural frequencies can be determined 
in any order desired, again just as in Holzer’s method. 
It should be noted that influence coefficients are not re- 
quired for the calculation and, indeed, they may be 
easily obtained by proper use of this method. In ap- 
plication to airplane wings, a single calculation for one 
value of K suffices as a trial for both symmetric and 
with only the last matrix 


[Ay] X [By’] 


[Fy + i] x 
[Fv] X [By’] 
antisymmetric modes 
[Aete ona) changing. 
TORSIONAL CALCULATIONS 
The equivalent associated matrix for the torsion of a 
beam section is (Fig. 5): 


‘T'y +1 1 KJy Ty’ 
both ends of th 6” 7 = jomeee Tyo 
oth ends of the nNt+1 Cy Cy N 


This equation was given in a somewhat different 
form by Pipes.’ 
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CHANGE IN DYNAMIC CHARACTERISTICS ACROSS COMPLETE BEAM SECTIOM 
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FIGURE § 
CHANGE IN TORSIONAL DYNAMIC CHARACTERISTICS ACROSS COMPLETE SHAFT SECTION 
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TABLE 1 
Typical Matrix for Beam Section of Wing 




















1 0 0 0.000472K 0 0.001444K 
60 1 0 0.02832K 0 0.08664K 
0.2184 0.007282 1 0.0001031K 0 0.0003154K 
4.369 0.2184 60 1 + 0.002062K 0 0.0063088K 
0 0 0 0.001444K 1 0.39999K 
0 0 0 —0.0000127K —0.¢ 008803 1— 0. 0035211K 
TABLE 2 
“Collapsed” Matrix for Entire Semispan, Antisymmetric Case, a= 1.8 
Vo ; bd 16. 181646 0.0763676 : —0.4555592 Vio 
Mo bh ui 2,315.4434 11.15001 —84.95284 Mi 
do { 7 .3919922 0.0363667 —0.1892385 dio 
yo = are : 1,490.4798 7. 1418962 ‘ — 16.039875 x 10 
To prs —261.97785 —1.2081959 65.339157 Ti 
4 none 1.3158321 0.0059174 0. 5484001_ A10 


~ Columns | 1,2 2, >, and 5 are not required for ‘the frequency solution; hence, they need not be computed, 
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For K#/.8 For K-.6 
O~M #23/5. $434 +11. /5001- 84.952840,, (1) : ~V 22. 5597363 $. +. O7O2ZESS ~.997F253@,, (i) 
0 = y = (490.4798 $, > 7./4/8962 - 16.0398756,, (2) ® 6 2.976225 9 +.009492 -.07672030,, (a) 
0 = 6 =/.31S832!9, +.00S9IT+ +, 548400/6,, (3) : © Ts - 483. 910/35 B= 1.4.3 0/544 +712, 9831154, (3) 
From (2) end(3) $= -.0047843, 6,,:.000689/ From (2) end(3) $,.* 003/934, 6,, =~. 000/616 
*. From (1) M =. 0139693 . Fromll) V= -.0016/56 
FIGURE 6 FIGURE 7 
IMBALANCE OF END COND/TIONS IMBALANCE OF END CONDITIONS 
ANT/-SYMMETRIC Mvs Kk ; SYMMETRIC Vivek 








COUPLED BENDING-TORSION 


Finally, the complete matrix equation for the beam section for the coupled bending-torsion condition is: 





[ Vyta t 1 0 0 KWy 0 KSy | | Ve’ 
M'y + 3 by l 0 K Wrly 0 K Sub 
; ly? mt KWyly? KSyly? 
¢'y —— - 1 —— ( ? 
sil 2Fl, —— Ely DEL y , 2Ely On 
, = ly*® by? KW, vly* KS, vly* 
yy — — bw 14+— 0 xX | Ww 
aligthe, 6ETy 2EIy . t 6E Ty 6E “2 F 
T'y +1 0 0 0 pa 1 } to Ta’ 
, KSy ] KJy 
Oy +1 | 0 0 0 -— -—— 1-- a" 
ie , a Cy Cy Cy | = . | 
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In practice it has been found advantageous to divide 
both inertia and stiffness terms by 10® in order to ob- 
tain a more uniform range of numerical values. Table 
| shows a typical matrix for a beam section of an air- 
plane wing carrying heavy concentrated weights. 
Table 2 shows the ‘‘collapsed’’ matrix for the entire 
wing for the antisymmetric case for a value of K equal 
to the Ist mode natural frequency. Fig. 6 is a plot of 
the imbalance of Jeftsiae versus K. A_ typical 
calculation of Mie siae is shown below the figure. 
Fig. 7 shows similar results for the symmetric case. 

METHODS FOR ABBREVIATION OF CALCULATIONS 

As the right side of the vibrating wing is always con- 
sidered free, Wright side» Might sider ANd T” sight side 
are always identically equal to zero. This means that 
columns 1, 2, and 5 need not be computed since their 
elements are always multiplied by zero. 

In practice it will be found that it is preferable not 
to combine the elastic and mass matrices before multi- 
plying by the trial value of K. This permits each of 
the premultiplying matrices to have a maximum num- 
ber of terms either zero or unity. 

The computation of a trial value of K actually re- 
quires considerably fewer operations than are required 
by a single iteration of the equivalent matrix by the 
iteration process. The actual time required by ordi- 
nary calculating machine, however, is about the same. 
This is-caused by the numbers having a varying order 
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of magnitude. When one considers that each trial 
value of K suffices for both symmetric and antisymme- 
tric cases, that no elimination matrices need be found, 
that relatively few trials need be made for any fre- 
quency, that influence coefficients need not be deter- 
mined, and that the modes may be determined in any 
desired order, the merits of the methods of this paper 
become obvious. 

The computation of an uncoupled bending mode of 
a ten-mass system for one trial value of K takes some 
40 min. by ordinary calculating machine. This means 
one value of K for symmetric and one for antisymme- 
tric. 

In determining the coupled modes, work was per- 
formed on the I.B.M. Here, a single trial value takes 
about 3 hours (ten masses). This involves, of course, 
no engineering time. By calculating machine approxi- 
mately 5'/. hours were required. 
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Letters to the Editor 


(Continued from page 566) 


In discussing the inclusion of the kinetic energy of the fuel, Mr. 
Szczeniowski asks: ‘‘Why not include this energy (the internal 
energy or molecular kinetic energy) also?’”’ This was done, 
although in a form that, perhaps, was not recognized by him. 
In any combustion calculation and, naturally, in any jet device 
involving combustion, not only the heat of combustion is con- 
sidered but also the internal energy of the fuel and air and the 
flow energy of the fuel and air, as well as the energy of molecular 
dissociations. All of these factors are included in the determi- 
nation of the jet velocity, which is, of course, a function of the 
thermal efficiency. The intra-atomic energy for any atomic- 
powered propulsion unit will also be a function of the thermal 
efficiency. 

Since confusion okviously exists as to the instantaneous pro- 
pulsive efficiencies of the various systems of propulsion by reac- 
tion, it seems desirable to summarize the characteristic features 
This summary is given in Table 1, page 566. 

M. Z. KRZYWOBLOCKI 

R. W. McCoy 
Associate Professors 
Aeronautical Engineering 
University of Illinois 


of each. 





Dear Sir: 

I wish to add the following remarks to the preceding letter by 
Professor McCloy and myself with reference to Mr. Szczeniow- 
ski’s Letter to the Editor (JouRNAL OF THE AERONAUTICAL 


SCIENCES, p. 310, May, 1947). In our previous letter was writ 
ten: ‘Assume that the mass of air has a negligible velocity with 
respect to the propelled unit as it enters the combustion cham 
ber.”’ It is clear that the kinetic energy (wa/2g) V? must be im- 
parted to the air. This quantity of air initially at rest with re- 
spect to the ground, must, when brought to rest with respect to 
the combustion chamber, be instantaneously moving with the 
velocity, V, with respect to the ground. This kinetic energy, 
imparted to the air by the propelled unit, is utilized as work in 
forcing that quantity of air into the combustion chamber against 
the combustion chamber pressure which is caused by ‘“‘choking” 
in the exit to the combustion chamber. After energy addition 
due to combustion, the air, in the form of products of combus- 
tion, appears in the exhaust with a velocity greater than V. As 
the basic period of time, 1 sec. was selected. In reality, the 
process is continuous, and unit time was used only for conveni- 
ence. Hence, it is obvious that all the calculations were per- 
formed with respect to one and the same reference system—i.e., 
the ground. 

In emphasizing the energy that must be imparted to the air as 
it is accelerated from rest with respect to the ground to the veloc- 
ity V as it enters the combustion chamber, the derivation was 
repeated with the propelled unit at rest and the entering air 


moving at a velocity V. 


M. Z. KRZYWOBLOCKI 
Associate Professor 
Aeronautical Engineering 
University of Illinois 
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Stability of Rotor Blade Flapping Motion 
When the Hinges Are Tilted. Generalization of 
the “Rectangular Ripple’ Method of Solution 


GABRIEL HORVAYi? anp S. W. YUAN? 


McDonnell Aircraft Corporation 


ABSTRACT 


The purpose of this paper is twofold: (1) to extend the well- 
known “rectangular ripple’? method of solution (Meissner, Van- 
derPol, Strutt) of certain differential equations to damped sys- 
tems with both first and second harmonic fluctuations in the coef- 
ficients; and (2) to apply the method to the determination of the 
flapping transient of a helicopter rotor blade. 

Several writers have been preoccupied with the nature of the 
transient solution of the flapping equation, since this has a bearing 
on the rapidity of control response. Glauert and Shone solved 
the equation by neglecting the inertia term and obtained the re- 
sult that the motion is unstable. J. A. J. Bennett (Aircraft 
Engineering, May, 1940; he also gives an account of Glauert’s 
work) solved the problem for the special case 6; = 0, » = 1.0, 
and m = 1.5 and found that the solution is stable. 

In the present paper a more general approach is adopted. It 
is shown that the stability problem of blade motion reduces 
essentially to the determination of two functions: (a) the natural 
frequency © of blade flapping, and (b) the apparent damping coef- 
ficient #app. of blade flapping. Instability cannot occur unless the 
blade frequency is an integral multiple of half the rotor frequency 
wand unless Mapp. is less than zero. It is found that the flapping 
frequency in forward flight is prevalently 0 for 6; ~ —15°, 
'/ym for 6; ~ 0°, and w for 6; ~ 23°, Yet, because of the large 
aerodynamic damping that a blade experiences during flapping, 
instability cannot set in except at unattainably large velocities of 
advance or at negative 6; of about 15° and over. 


NOTATION 

w = rotor frequency 

t = time 

¥ = ot = azimuth (measured from rearmost posi- 
tion of the blade) 

7 = normal flapping angle (measured from the 
hinge axis) 

8 = flapping angle (measured from the un- 
deflected position) 

8 = dB/dt 

8° = dB/dy 

B* = conjugate complex of 8 

R(8) = real part of 8 

§(8) = imaginary part of 8 
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17, 1947, 
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Q(x) 


n = pcaR*/8I 
Napp. 
u 


Napp./n = 1 
Napp./N <0 

1 — Mapp./n 
F(8) = E(wt) 
C(t) 


P*(t) 
E(t) 


90° — 43 


ry 
Qa 
U, Up, Ur 


C., Pe, G, Po, etc. 
Pp 

pi’, p:?, ps’, pi? 
Ba, Bo, etc. 


Q 
N 


e(y) 

1) 
‘i, te, Se’, Ba 
Ai, B,, As, Bo, ete. 
0, Pc, Ox, Ore, O2s 
£0, £1, hy, Go, Gi, Hy 
Qk, by 

Oy) ~ 


natural frequency of flapping (denoted by 
o,w in reference 4) 

fractional part of Q/w(—1'/2< 7 < +1/2) 

flapping frequency at » = 0 

advance ratio 

inflow ratio 

manually set blade pitch angle 

actual blade angle 

angle of attack 

air density 

chord length 

rotor radius 

moment of inertia of blade with respect 
to the flapping hinge 

aerodynamic parameter 

apparent damping coefficient 

1/3(m — Mapp.) (denoted by a; in reference 
4) 

criterion for ‘“‘complete’’ stability 

criterion for instability 

degree of destabilization 

equation of flapping motion 

coefficient of 8 in $(8), defined by Eq. 
(1b) 

coefficient of 8 in $(8), defined by Eq. 
(1c) 

aerodynamic excitation, defined by Eq. 
(le) 

angle that the flapping hinge forms with 
the: blade axis 

function defined by Eq. (2) or (44) 

lift force 

centrifugal force 

mass of blade element 

distance of dm from the flapping hinge 

distance of dm from the hinge axis 

arc of travel of dm during flapping (Fig. 1) 

air velocity and its perpendicular and 
tangential components 

certain average values of C(t), P(t) 

C/w, P/w 

Pa? von ¢*, pr? a Co, pe? oan ce. pa® _— Ca? 

representatives of 6 in certain azimuth 
intervals 

factor of amplification 

coefficient in the consistency criterion, Eq. 
(20) (denoted by 1 — 2g? in reference 4) 

exponential factor of 8 

oscillatory factor of 8 

quantities defined by Eqs. (20a) and (34) 

coefficients of Ba, 8», etc. 

coefficients in Eq. (3) 

coefficients in Eq. (48) 

coefficients in Eq. (51) 

“except for a normalization factor, ®(y) 
is equal to”’ (Fig. 3) 
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INTRODUCTION 


, I ‘WE DIFFERENTIAL EQUATION of flapping motion of 
- a rotor blade is, as shown in the next section, of the 


form 

5(8) =B + 2C()8 + P28 = E(ot) (1) 
where 
E(wt) = nw? {4/s\ + (1 + w2)d8 + 


(2h + °/3u) sin wt — wd cos Qut | (la) 


is the aerodynamic forcing function and C(t) and P?(f) 
are periodic coefficients, 


2C(t) = nw(1 + 4/3u sin wf) (1b) 
P?(t) = w?[1 + mu(4/3 cos wt + w sin 2wt) + 
n tan 63(1 + w? + 8/3u sin wt — pw? cos 2wt)] (lc) 


The substitution 
BY) = vy) X exp. (—'/2 mb + ?/snp cos y) (2) 
transforms Eq. (1) into 
d*v/dy? + [00 + 0.cosy + 0; sin Y + 2, cos 2y + 

Bo sin Qy le oe e /ony — */3np cos yp E(wt) (3) 
where 
0 = 1 — !/qn? — */gu?2n? + (1 + w?)m tan 43 
6, = —?/sn2u + 8/snu tan 43 


2/on*u? — ny? tan 63, Oo, = np? 


6, = */snu, 
6. = 


The value of w?% at uv = 0, 
Q? = w?A(0) = w(1 — '/qn? + n tan 63) (4) 


is the squared natural frequency of flapping at » = 0. 

There are two questions that arise in connection with 
Eq. (1): (1) What is the forced response 8) of the 
blade to the excitation E(wt)? In particular, is reson- 
ance excitation possible? (2) What are the two natural 
modes 6; and f» of the blade—i.e., what are the two 
solutions of the homogeneous equation 


F(8) = B + 2C(t)B + P*(t)B = 0 (5) 


In particular, is instability possible? The question 
of resonance arises from the observation that, for suit- 
able values of 6; and yu, the blade frequency Q can as- 
sume the values kw(k = 0, 1, 2, ...), where w is the 
frequency of the exciting function E(wt). The question 
of instability arises from the observation that, for suit- 
able values of 6; and y, the blade frequency 2 can assume 
the values 1/,kw(k = 0, 1, 2, ...) where w ts the fre- 
quency of the periodically varying coefficients C(t), P*(t). 
The first set of 2 values makes the flapping motion 
susceptible to resonance excitation; the second set 
makes the motion susceptible to instability. 

The two questions lead to two intimately connected 
problems of: (a) determining the natural frequency 
of blade flapping 
(6a) 
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as function of the advance ratio 4 and the parameters 
n and 63; and (b) determining the apparent damping 


coefficient of blade flapping 


(6b) 


nN = Napp. (qu) 


app. 


as function of » and the parameters » and 6;. The 
terminology ‘‘apparent damping coefficient”’ is justified 
by the observation that, while the equation 


(d?B/dy?) + (ndB/dy) + (1 + n tan 6:)8 = 0 (7) 


to which ¥(8) = 0 reduces when pu = 0, has the solu- 
tions 

8, = exp. (—1/snwt) cos Qot \ (2) 

, J J ) 

8B. = exp. (—!/2nat) sin Oot § ‘ 


the more general equation 


+(8) = O, pw #0 


has, as will be shown later, solutions of the form 


By = e,(t)P,(24) 
Bo = es(t)Po(Qt) § 


where &; and &, are oscillatory functions of frequency 


Q, and 


e,(t) = exp. (—!/2Mapp.,wt), 


e2(t) = exp. (-! ‘oMapp.,wt) Ya) 


Napp., = N + 2u, 


u(u, n, 63) > O 


Napp, = nN — 2Uu, 
(9b) 


Evidently, mapp., is the critical damping coefficient. 
It will be written without the subscript “1.” When 
Napp. = nN, the motion f; (and, hence, also the motion 4, 
which is a linear combination of 8, and 82) is completely 
stable; when m > Mapp, > 0, the motion (still stable) is 
partially destabilized; and when Mapp, < 0, the motion is 
unstable. 


DIFFERENTIAL EQUATION OF BLADE FLAPPING 


Eq. (1) will now be derived on the basis of the follow- 
ing assumptions: The blades are attached to the rotor 
shaft by means of flapping hinges only (the freedom of 
oscillations of the blades in the rotor plane about the 
lag hinges is thus disregarded), and the distance of the 
hinges from the rotor axis is negligible. The hinge axes 
are perpendicular to the rotor shaft but may form the 
arbitrary angle 90° — 63 with the blade span axes. The 
blades are untwisted and are rigid both in bending and 
in twist. The chord and the lift curve slope are constant 
along the blades. The uniform motion of the helicopter 
is not affected by the oscillations of the blades. Of the 
aerodynamic forces that act on the blades, only the 
velocity square lift forces are taken into account. Air 
velocities in the direction of the blade span are neg- 


lected. The field of induced velocities is assumed as 
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STABILITY OF ROTOR 
constant over the rotor disc, and interference effects 
between the various blade elements are disregarded. 
All angles involved (but Y and 63) are small; in par- 
ticular the angle of attack ts always less than the stalling 
angle. 

For purposes of the present analysis, only constant 
blade angle settings are considered. When the blade 
flaps, at angle 


B = y cos 63 (10a) 


with respect to the rotor plane (y is the flapping angle 
with respect to the hinge axis), the blade angle reduces 
to 


ov = 9 — Ad (lla) 


where, cf. Fig. 1, 


Ad = a/r cot 63 = y sin 63 = B tan 6; (11b) 


The equilibrium of the air, deadweight, centrifugal 
and inertial moments with respect to the hinge axis 
gives 
Seer, (dL — gdm — ydC cos 63) — 

ee 
JSr=o ¥r,"dm = 0 (12a) 
Here, 


r, = rcos 63 (10b) 


is the distance of the mass element dm from the hinge 
axis (r is the distance measured along the blade span 


axis), 
dC = rw*dm (12b) 
is the centrifugal force acting on dm, and 
dL = '/.pca(x + 8’) U°dr, (12c) 


is the air force acting on dm. 
Noting that! 


Ur = rw + woR sin y | 
Up = wR — 1B — pwwRB cos yy? (13) 
(x + 8)U2& UpUp + 8'Ur? { 


applying Eqs. (10) and (11), performing the integra- 
tions indicated in Eq. (12a), and neglecting the dead- 
weight moment, one obtains the flapping equation in 
the form of Eq. (1), Q.E.D. 

There are many ways in which a tilted flapping hinge 
can be achieved kinematically. Recommended by its 
simplicity is the solution indicated in Fig. 2, which 
shows a control horn displaced at angle 6; from the 
normal flapping hinge. 


The pair of differential equations so obtained 


A + 26.8’ + Pa*Ba — 0 
By” + 2cyBy’ + Po Bo = 0 
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SOLUTION OF THE HOMOGENEOUS FLAPPING EQUATION 


Twin-Ripple Solution®* 


The aim is to find the solutions 8; and 8: of the homo- 
geneous Eq. (5). Instead of considering the exact prob- 
lem, 5(8) = O, with the functions C(t) and P?(t) as 
given by Eqs. (lb—lc), one can consider the simpler, 
approximate problem, II,,, which is obtained when one 
replaces C(t) and P?(t) by their time averages C, and 
P,” in the azimuth range Yq to Yo + m and by their time 
averages C, and P,” in the azimuth range Yo + 7 to 
Yo + 2x. For instance, in the approximation II_,/2, for 
63 = 0, 


ZL. = 2C, = nw ( 
P,? = w*[1 + (Sny/3z) | sales 
fe = w*[1 — (Sup /3m) i) 


and in the approximation IIo, for 6; = 0, 


2C, = nw[1 + (Su/3m) ]) 
2C, = nw[1 — (8u/3x)]> (15) 
P2= Pf =~ \ 
foro < wt C y+ t (16) 
for Wo + wr < wt < wo + 2rf , 
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where 
8’ = dB/dy, pb = P/», 


have the solution 


B= ‘ Ba = eA; cos pw + B, sin piv) 
By = e “(As cos poy + By sin poy) 
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c = C/o, 


1947 


Here, A;, B,, Ao, By are determined from the conditions 


Bo(Yo + m) = BalYo + =), 
By(Yo + 27) = QBa(Yo), 


pi? = pa? — Co", = ra? = Py? — Oy? (17) 
Yo g v g Yo +f 

votridvSwt+2n (18) 

By’(Yo + 7) = Ba'(Yo + 7) (19a-19b) 


By’(Yo + 2m) = QBa’ (Yo) (19c-19d) 


which express continuity of function and derivative at the junction azimuth y+ 7, and a Q-fold increase in the 
otherwise periodic 8 and §’ in the course of a complete cycle 27. 


Substitution of the expressions (18) into the homo- 
geneous system of Eqs. (19) indicates that a solution 8 
of the desired form is indeed possible provided the factor 
of amplification per cycle, Q(27), is a root of the deter- 
minantal equation 

(Q/€a€)” — 2N(Q/e.@) + 1 = 0 (20) 
where 
N = cos rp; cos rp2 — 
pr? + P:° + (c. — &)* 
2pipe 


sin rp, sin rp, (21) 





and 








— »~cat —chr 
é 


ins, q = , €6 =e "* (20a) 


One thus obtains the result that for 
Q = (N + VN? — le" (22) 


the Eqs. (19) are consistent and can be solved for the 
coefficient ratios A,/B:, B,;/B2, A,/B,. (Evidently one 
of the coefficients, or what amounts to the same thing— 
the normalization factor of 6—can be selected arbi- 
tarily.) Two sets of coefficients are obtained, accord- 
ing to whether the plus or the minus sign is used with 
the square root. They define the more stable solution 
8B. and the less stable solution f;. 


It will be convenient to define functions u and r of u, n, 53 by the relations 


Gi) &™=N+ SN? — lforN>1 ) 
(Gi) e"* = —N+ WN? — 1 for N < —1 ( (23) 
(iii) 2x7 = arctan x/1 — N?/Nfor -—1< N <1 
Then 
(i) an = ” le cata eal : QO, a ela /m—u)en for N > 1 
Gi) Dy = gf eter O, «a flee" tae N < 1 (24) 


=x ef—'/m— i)" for —1 € NQ 1 


~) 


, 
Git) Q, = ef vertinee, Q 





(—Mapp,¥) and of the oscillatory function 


are the amplification factors per revolution of 8; and exp. 
@, is periodic when 7 is a rational 


Be. One notes that u > 0 in cases (i) and (ii) and that @,(~) = Bi(W)/a(y). 
u = Qin case (iii); furthermore, 7 = 0 in case (i), 7 = number. 

+1/, in case (ii), and —'/2 £ + < 1/2 in case (iii). In particular, ®, is integral periodic for r = 0 and half 
Thus, the expressions for Q indicate that there is both periodic for 7 = 1/2. (It is convenient to use the ex- 
pression; ©, is “purely periodic’ for 7 = 0 or 1/:.) 


real amplification (or reduction) 
Corresponding results also hold for 6. 


(25a) 
The above analysis yields an explicit expression for 


exp. (—Mapp., 1) 





involved, as well as imaginary amplification (i.e., a 
phase shift), 
(25b) 


exp. (27772) 


The forms (9) are now established: ; is in effect 


the product of the exponential function e(y) = 


the apparent damping coefficient m,pp, in terms of 
N(u, n, 63). The fractional part + of the companion 
function, the frequency ratio, Q/w, is also obtained 
directly, but the integral part becomes evident only 
after determination of the function pair 8,, 9p. 
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Twin-Ripple Solution for Double Frequency Fluctuations 


For 6; = 0, uw > 1, the uw terms of Eqs. (5) can be 
neglected in first approximation, and one obtains 


B+ nwB + w? (1 + np? sin 2wt)B = 0 (26) 


This equation no longer represents flapping motion 
(helicopter operation at 4 > 1 is meaningless), but the 
mathematical problem of solving Eq. (26) is still of 
interest. 

The method of the preceding section can be readily 
adapted to this case as follows: Consider the approxi- 
mation IIy’, where the subscript ‘‘0’’ refers, as before, 
to the value of Yo from which the first ripple is counted 
and where the superscript ‘‘2’’ points to the fact that 
two complete fluctuations take place in the 27 interval. 
Denoting the averages of the 8 and 8B coefficients by 
2C,, P,” in the interval 0 to '/27 and by 2C,, P,? in the 
interval '/27 to , the solution 8 of Eq. (26) can again 
be represented by functions 8, and £8, for the respective 
intervals, as in Eqs. (18); but now 


'/49) = B,'(!/er) (27a-27b) 
B,’(r) = QB,'(0) (27c-27d) 


By(!/2r) = Ba(1/em), 
B(r) = Q8,(0), 








Four-Ripple Solution 


BLADE 
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are the boundary conditions from which the coefficients 
A,, B,, Az, B, must be determined. The consistency 
determinant of Eqs. (27) assumes the form 


(e'""Q)? — 2N(e’"""Q) + 1 = (28) 
where 
N = cos !/erp; cos '/erp. — 
pi* + pot + (Ca — &)? 
2pipe 


The amplification factor per cycle of fluctuation (i.e., 
per azimuth interval y = 7) is therefore 


, 


sin '/orp; sin 1/emp2 (28’) 


= (N = ~/N? — I)e'"" (29a) 
This can be written in the form 
Q = exp.(—!/on = u + ir) (29b) 


Correspondingly 6, and 6, factor into ¢; X &, and 
é2 X Po, wheree, ,(y) = exp. (—}, ‘Mapp., _y) and @,, X 
(QY/w) is integral periodic for 7 = 0 and half periodic 
for 7 = l. 

The extension of the method to higher frequency cases 


is obvious. 


Subdivision of the basic y interval of 0 to 2x into four equal sections (approximation IVo) and replacement of 


5(8) = O by the equations 


Ba” + 2¢aBa’ + Pa’ Ba = 0 0 g y g a/2 
By" + 2crBy’ + pr?By = 0 m/22¥Qg4 (30) 
a + 2c.B,’ + P-*Be = 0 us g y < 32/2 : 
Ba” + 2ceBa’ + pa’ Ba = O 3/2 < v $20 
obtained on averaging 2c(y), p*(W) over the respective intervals, leads to the solution 
Ba = eA, cos pw + B, sin pw) 0 g v < n/2 
3 B, = e “(As cos poy + Be sin po) m/2Rvg4 (34) 
~ |B, = e~“*(A; cos pw + Bs sin psp) tiv < 3x/2 , 
Ba = e “(Aq cos pa) + By sin py) 39/2 Sy < Qe 
subject to the boundary conditions 
By(4/2) = Ba(w/2), By’(r/2) = Ba'(m/2), Br) = Bm), B.'(m) = By'(r) — (32a-32d) 
Ba(3x/2) = B,(34/2), Bq’ (34/2) = B,'(34/2), Ba(2r) = Q6,(0), Ba’(2r) = QB.’(0) (32e-32h) 


Approximation IV» is an important generalization of the twin-ripple method. 


It permits simultaneous evalua- 


tion of both first and second harmonic fluctuations in the damping and spring constants of the system. 


Eqs. (32), written in full, are: 





€C;! €a51! — 6C —_ €So! 0 0 
€a(— pi S;! = CC") €a(piCy} —_ CaS) stan + °C) €o(— p2oCo} + CoS2") 0 
0 0 2? €7.S2? —€,2C;? — €c7S3? 


€c?(— ne + ceSs*) 


0 0 ops 2 — ¢4C>?) €5?(p2C2? — cyS2*) eo*(DeSs 2+ c.C;*) 
0 0 0 2 C33 3 S3 
0 0 0 0 e.(—paSs 3 — ¢.C;3) MAG — ¢S33) 
Q 0 0 0 0 
—Qca Op: 0 0 0 0 
0 0 | | Ai 
0 0 | By 
0 0 | | As 
0 0 1X] By i= 0 (33) 
—€8C,3 — €43S,3 ( | As 
€a*(ps SGtic aC;°) €a°( — psC,3 + caSy3 ) | | B, 
—egtC,! —eqiSy | | Ag 
€a*(paSgt + cay’) — €a*(— pag + caSe*) | |B, 
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Here 


9 


er — ea die ir = a,,0; ¢ a) 
C’ = cos jp,x/2, SZ = sin jp,r/2 


t (34a) 


The determinantal equation of (33) is again 


(e""Q)? — 2N(e""Q) +1 =0 (20) 
where, using the notation 
C, = C;, Dy, = Sy'/pr (34b) 
one can write 
2N = 24,0.03C, — D,D2C3C;[p1" + p2? +(c, cat Cy)” ] ial 
D,C2D3C;[p1? + ps? + (Ca — €-)?] — 
C:D2D3C[po? + ps? + (@ — ¢-)?] — 
D,C2C3D4[pr? + ps? + (Cg — €a)?] — 
CiD2C3D4[po” + pa? + (Go — Ca)*] — 
CiC2D3Da[ps* + ps? + (¢e — Ca)?] + 
D,D2D3Cs{(pr2 + e¢e) (Cp — Cc) + 
(pe? + Gila) Co 7 Ca) + (p3” + Caln) (Ca ial Co) | + 


D,D2C3D4[(pa® + Cale) (Ca — Cro) + 
(Pi? + CoCa) (Cy — Ca) + (D2? + Cala) (Ca — Ca) |] + 
D,C,D3sD4[(ps? + Cala)(Ca — Ca) + 
(pa? + Cale) (Ca — Cc) + (pr? + Cc6a) (Ce — Ca) ] + 
C:\D2D3D4[(p2? + CeCa) (Co — Ca) + 
(ps? + Caly) (Ca — Co) + (pa? + nee) (Cy — co] + 
D,D2D3Di[pir-ps + po" pa? — Ca7C,? — Cy7€a” + 
Walrela + (pr? + Ca”) (C.? — Cely — Cela + Cola) + 
(P2® + Cy”) (Ca? — Cala — Cale + Cale) + 
(bs? + €.”)(Ca” — Caly — Cala + Cola) + 

(pa? + Ca?) (Cy? — CyCq — CoCo + CaCe)] (35) 


For the normalization B, = 1, one obtains 


The factors e;(y) =exp. (—!/2Mapp.¥) =exp. (— 0.5893 X 
vy) and P,(Y) = Bi(v)/er(y) of 6, are drawn in Fig. 3a,* 
the first normalized to 1 at Y = O and the second nor- 
malized to 1 at its maximum. The determination of 
the solution 82 is similar and is also shown in the figure. 
In Figs. 3b and 3c* the same curves are shown as 
obtained by the methods IIp and [Vo. It is seen that 
the natural frequency. of flapping, 2, is one-half of the 
rotor frequency w. 

The close agreement of the m,pp, values in the various 
approximations is at first rather surprising. Offhand, 
one would expect that the fluctuations in the damping 
constant (evaluated by approximation II») and the 
fluctuations in the spring constant (evaluated by ap- 
proximation II_,,/2) are additive in their effect and 
should result in a diminished n,,,, when both effects are 


* In Figs. 3a and 3c the labeling of the curves e and e; should 


be interchanged. 


8 Ba = e °**¥(0.4798 cos 0.8825y — 0.5776 sin 0.8825y) 0 
' 1 By = e~°**(0.2583 cosh 0.4731y — sinh 0.4731y) ~ 


SCIENCES—OCTOBER, 1947 
The functions « and 7 are defined as before by Eqs, 
(23), and the interpretation of the results also remains 


unchanged. 


Application to the Flapping Problem 


It will be profitable to solve the flapping equation by 
methods II_,/2, Ilo, 1Vo, and IIo” for the representative 
parameters m = 1.7, uw = 0.34738, and 6; = 0. This 
will permit a better appraisal of the relative merits of the 
various approximations. 

By Eqs. (14), (17), (21), and (22), one finds that in 
the II_,,/2 approximation 


26, = 2% = 1.7, pr = 0.8825 
p, = 0.4731i, N = —2.6694 (36) 
Q,(27) = —0.02467 = — exp. (— 1.17867) 

Q.(27) = —0.0009322 = — exp. (— 2.22127) 


and therefore the Eqs. (19) for the coefficients A, By, 
A», Bz of B, assume the form 
+ 0.9830B, — ) 
1.2890A» 
— 0.6736B, + 
0.7109A>2 
0.01721A, — 0.09218B, + 
0.08557A2 + 0.08361B7 = 0 | 
0.06672A, — 0.09354B, + | 
0.033184, + 0.03059B,i = 0 } 


0.1836A, 


— 0.8134Bt 0 


—1.0235A, 


| 
+ 0.0816Bi = 0 ean 
r (37 


No a 
3 
—_—— 
_ 
y ~ 
2 


Z2y2 
Z2w2 


taken into account, as is done in approximation IV. A 
closer investigation, cf. Eq. (45), confirms the observa- 
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tion that the two effects partially offset each other. 
On the other hand, the functions ®, and ®, obtained by 
the three methods differ greatly, and this indicates that 
the twin-ripple approximations IIp and II_,,/2 are 
not suitable to represent the wave form of blade 
flapping. 

Solution of Eq. (5) (for w = 0.34738) by method IT)? 
is not expected to yield functions 6 that resemble the 
correct ones, because II, is applicable to Eq. (5) only 


when uw» > 1. But it is useful to include also a solution 
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behavior of non-purely periodic functions. 
By IIo, one has 


by IIo? among the examples, because it illustrates the 

26, = 24 = 1.7, pi = 0.63883, 

p2 = 0.38327, N = —0.098405 20) 

O.(r) = e °**(—0.098405 + 0.995157) = | at 
exp. (—0.85a + 95.65°7) = Qo*] 


.= 
and therefore the Eqs. (27) for the coefficients A, B,, 
As, By assume the form 


e °°" 0.537384, — 0.84334B, + 0.824184, + 0.56633B, )= 0 
e~ °#5*( _9.99552A, — 0.37355B, + 0.91761A2 + 0.16549B:) = 0 jm 
e 8 9 S704, + 0.358542 + 0.93351B2) = 0 | | 
e °**"(_0.85e°"0A, + 0.63883e°*"0B, + 0.66255A2 + 0.65607B:) = 0 
With the normalization A; = 1, one obtains the solution 
{ 8, = e 9 ¥ [cos 0.63883y + (0.2119 + 0.82277) sin 0.63883y] OS < Yor (41a) 
b= , 
(By = e~ (1.2788 + 0.155452) cos 0.38327p + (— 0.5967 + 1.00634) sin 0.383279] '/ar Sy <x (41b) 


for Q1, and By = 8,* for Q:. 


But the natural modes of oscillations are necessarily real. One obtains the two real transients 8; and 3 by 


taking 
R(B,) OS ¥ << #/2 
R(By) 4/2 g ¥ NE 
R(QB.) - Jr iv < 3x/2 
B. = for \,_, ; 
R(QBy) 3m/2 Qy < Qe 
&(Q78,) 2e Lv < 52/2 


(It is really immaterial which function is called 8; and 
which is called 62, because the extinction factor exp. 
(—"'/oMapp. ¥) = exp. (—O0.S5y) is the same for both func- 
tions. Therefore the convention of calling the less 
stable solution @; is irrelevant in the present case.) The 
functions ®, and @, are plotted in Fig. 4 for the range 
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5 (Bq) O<y¥C 2/2 
I(By) m/2RWKQr 
. 5(QB,) tiv < 32/2 
8, = , 2 
= \5(0m) “ \3x/2 < ¥ < 2m 
5(Q28,) Qe < ¥ < 5e/2 
0° < ¥ < 1,800°. They were calculated in 45° steps. 


It is seen that ®, and @, are oscillatory but not periodic. 
The pattern of the wave shapes @), 2 varies indefinitely, 
without repetition, as y varies from —© to +. In 
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particular, it can be noted that the @, intersections with 
the abscissa axis occur at about 0°, 339°, 680°, 1,020°, 
1,361°, 1,700°; the ®, maxima, at about 196°, 846°, 
1,514°; the ®, minima, at about 524°, 1,179°. The 
arguments for ® are similar. The inconsistent varia- 
tion in the cycle length gives rise to the impression that 
®, and @®, are not associated with a fixed frequency. 
Nevertheless, the functions @; and @2 have a definite 
frequency; from the expression for Q(7) one sees that 
this is 2 = 95.65w/180 = 0.5314w. 

As a concluding example, two cycles of the unstable 
V/3, 63 = 0, » = 2 are shown in Fig. 
In this case, 2 = w. 


transient, m = 
5 as calculated by method IV. 


Comparison of the Approximations 


On the basis of what has preceded, one can designate 
the most suitable ripple method of solution of a differen- 
tial equation with fluctuating spring and damping coef- 
ficients. 

Method II can evaluate the effects of first harmonic 
sine fluctuations only; method II_,/2, of first harmonic 
cosine fluctuations only. They give zero averages for 
all other harmonic fluctuations. Thus, IIo is used when 
C(t) and P(t) contain only sin y terms; II_,/2 is 
used when C(t) and P*(t) contain only cos y terms. 
Note that nothing is gained by using IV» in the above 
cases. 

Recalling that an equation of the form 


B” + (go + g. cosy + h sin y)p’ + 
(Go + G, cosy + Hi, sin y)B = 


can be transformed by means of the substitution 


0 (43) 
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B = v()[exp. (—1/ego — 1/2g1 sin Y + 
1/,h, cos ¥)] (44) 

into 
0” + [(Go — */ego? — */eg:? — */shy?) + 

(G, = 1/ohy sie 1 /sg081) cos y + 

(I, + 1/984 ma 1/o9oh;) sin y a 

1/3(hy? — gi?) cos 2W — 4/agiky sin 2¥]v = O (45) 


one readily sees that the more general Eq. (43) can 
also be solved by a twin-ripple approximation, provided 
that g," and h,” are small. In this case one can apply 
II,, to the equation 


vo” + [0 + V02 + 6,2 sin (Y — ¥o)]v = 0 (46) 


which is obtained from Eq. (45) by neglecting the second 
terms and combining the first harmonic 
sine and cosine terms. This method of solution 
can be called the II,; approximation. The bound- 
ary conditions for 8, Eqs. (19), reduce to the condi- 
tions for v 


harmonic 


Uy (Wo +r) = Va(Wo + 1) 
vy’ (Wo + 1) = V_' (Wo a 1) (47) 
*or(o + 2r) = Qua (Yo) 


ey,’ (Wo + 27) = Qu,’ (¥ Yo) 


and again lead to 
e"oO=N+VN- 1 
where JN is given by Eq. (21) in terms of 


C= % = 0 
pr? = 6 + (2/n)V 02 + 0,2 (48) 
p2* = % — (2/m)V 0.2 + 0.2 


When g,? and h,” are not small or when there are also 
second harmonic fluctuations in C(t) and P?(¢), then a 
four-ripple approximation must be used. Where the 
amplitudes of the second harmonic fluctuations greatly 
outweigh the first harmonic fluctuations, the latter 
terms can be neglected, and the approximation II? is 
applicable. 


Comparison with the Exact Solution 


‘Since both the twin-ripple and the four-ripple ap- 
proximations lead to the same form of the consistency 
criterion, Eq. (20) (only the value of N is different in 
the two cases), the question arises whether one may not 
obtain the correct expression for the exact N by a 
suitable passage to the limit, using infinitely many 
ripples. 

Indeed, an exact expression for N is readily obtained, 
and its expression is well known: 


Nexact = 1 — 2D sin’ rr/9% (49) 


where ® is the infinite determinant 
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= : : ‘ : < | 
y} (44) |. 1 0*G, 6*G, 0 O , | | 
| 4G, 1 6:*G; 62*G,; O seal ee SS APPRox mation ZT, , —_ 
D=|. Go AGo 1 A:*Go O2*Go .| (50) | pre 
|. O OG, 4G, 1 A*G, . | 
= ae G2 AG, 1 . 7 _——— 
0 (45) and I-M, ap /mn | 
43) can 0, = 0. + 16,, O2 = Br. + i025 \ sini | 
brovided G, = (0 — k?)— f (50’) 2 
n apply ' ’ 
: However, the procedure used in obtaining Eq. (49) is a | 
direct one and does not involve a passage to the limit, —— - 
0 (46) using ininitely many ripples. 
The usefulness of Eq. (49) was limited heretofore by 
> second = O 2 


the difficulty attendant to the evaluation of D. Re- 





armonic 5 : 
hati cently, this difficulty was removed by setting up tables 
solution P : 
ms for the coefficients of the first ten powers of uw in the 
O - ° ° . 
ae power series expansion of 9, covering a 4 range of 


-1.0to+0.9. Asaresult, the flapping problem is now 2 
accessible to exact solution in the range of the customary ; | 
values of w and n and for 63; between —15° and +15°. 
. Nevertheless, in the stability investigation of problems ~ 
(47) not covered by the tables, or, especially, in a study of 
cases where the power series expansion of D diverges 
(u 2 1), the rectangular ripple approximation remains 
an invaluable aid. 





























It will be profitable to compare Fig. 3c of the present 
paper with Fig. 2 of reference 4, which shows the exact 
flapping angles 8, and # for the condition n = 1.7, 
6; = 0, uw = 0.34738. It is seen that the four-ripple 
approximation of ®, and @, nearly coincides with the 





















































(48) exact ®; and @»., but the extinction factors e; and é» 

differ slightly : 

are also (Mapp./)4-rippie = 0.691 ’ (Stage, /% ennat = 0.638 

then a : : . : Fic. 6b. 

: Use of the R.M.S. value 1/./2 for the effective height 

ere the of the sin y loop in the formation of the equivalent 

greatly ripple, rather than the average value 2/7, brings the 

' latter approximate pp, into much closer agreement with the 

n II" ts correct Mapp. The writers now believe that R.M.S. 
values should be used in preference to average values in 
the formation of the parameters c,, p,’, etc. 

ple ap- STABILITY OF FLAPPING MOTION AT 63; = 0 

‘istency ; P : 

cout tx It was seen earlier that for a typical helicopter (nm = 

nay not 1.7) at the high speed of advance, u = 0.34738, the 

V by a degree of destabilization = 2u/n = 1 — Mapp,/n (51) 

— amounts to only 1 — (1.178/1.7) = 0.307 (by method 

tained, II_s/2, Fig. 3a). One therefore surmises that the flap- 
ping motion is very stable also at other values of u and 
n. This expectation is borne out by Fig. 6, which shows 

(49) the 2u/n values plotted vs. u for a variety of n. Areplot 














of Fig. 6 results in the stability diagram, Fig. 7. (Figs. ° 2 
6 and 7 were obtained by the II_,/2 approximation.) Fic. 7. 
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The shaded areas represent regions of complete stability; | 
the unshaded areas 0 and '/, indicate regions of partial bl 
destabilization. In the former, mapp./n = 1 (N is com- 
plex, of modulus 1), the rate of extinction of a transient 
is exp. (—!/any), and the natural frequency of oscilla- 
tion 2 is between 0 and !/.w; in the latter, 1 > mapp./n > on 
0, the rate of extinction of a transient is exp. (— '/2M%ap, at 
y), and Q is 0 (when N > 1) or !/2w (when NV <-— ex 
Since for most helicopters m is not far from +/ 3, the TI 
Fic. 8 operation of the craft at high forward velocities gen- Fis 
erally takes place in the Q/w = ‘/2 region. Jt is seen bs 
that for u < 0.5 the flapping is very stable. For instance, sta 
in the case of nm = 1/3, uw = 0.32, the transient extin- Fig 
guishes at the rate of 1 — exp. (—Mapp.7) = 1 — exp, V/ 
(—0.7\/3mr) = 97.8 per cent per revolution. ? 


Were the diagram extended to the right, one would use 
note that the lower boundary curve Mapp,/n = 1, which cat 
exits from the point n = 3, » = 0, reaches a lowest in t 
ordinate of about m = 1.34 at uw © 0.75 and then rises The 
again. Beyond this boundary curve there lies a region in I 







































































of complete stability, mapp./n = 1, in which the fre- T 
quency is between '/2w and w. To the right of this nati 
there follows the region Q/w = 1 in which mapp./n de- the 
creases from | to zero first and then to negative values. 
(Instability sets in. A particular case is shown in Fig. nb 
5.) Eventually, ,p,, rises again, a region of complete 
fA Fy stability is reached, and the process repeats indefinitely, This 
YL peer RonMaton Mag //)/ with @ increasing by !/2w across each region of complete | take 
re) a email ; 
icles stability. . . nha tion 
The approximate stability diagram, Fig. 7, shows oy 
good agreement with the exact stability diagram, Fig. 5 
of reference 4. 
i i : | > 
SY ) . ah, ms CONDITIONS FOR RESONANCE EXCITATION 
\ | ks . 
m NAC es % oe ~ | / The forced response §» of the blade to the excitation 
i a. ‘5 A | a | (wt) is of the form ae) 
18 . _ fm = 410 — 
: LS as T 
ioe LO | | “ = 
, , Bo = >> (a cos kwt + b, sin kot) (52 
an 8, = k=0 
V tank 2 VAP 7 | Resonance can exist only when the blade has a natural 
aii ogy. | ode 6; that is of the same f What are the re 
7 [APPROXIMATION [-% /| /> mode {; that is of the same form. at are 
{¢ rears ©, Yi Oe Pg. quirements for this condition? The transient 
0 22 “4 i +6 8 
Fic. 9b. Bip) = ex(t)P, (Lt) | 
G | 
can be of form (52) only when, first, e,(¢) = 1, or, what F | 
is equivalent, n,pp,/n = O—i.e., when the blade operates 
on the boundary of partial destabilization and instabil- 





ity—and, second, when the frequency Q of @, is # of 
an integral multiple of it. 

The lowest resonance of Eq. (1) therefore occurs 00 
the line m,pp./n = 0 of the Q/w = 1 region. On this 
line the steady-state amplifications are infinite. The 
u values for which large amplifications take place may 
spread well into the adjoining instability region on the 
right and into the stability region on the left. Since, 
Fic. 10. however, the #,pp,/n = 0 line runs far to the right of all 
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STABILITY OF BOTOR 
achievable advance ratios yw, resonance excitation of 
blade flapping is of no import in the present problem. 


ROLE OF THE 63; HINGE 


When the flapping hinge is cocked at a negative 6; 
angle < — arctan 1/n, then the excitation of flapping 
at hovering, exp. | {’,, exceeds the rate of damping, 
exp. (—'/2nwt). This follows from Eqs. (4) and (2). 
The region to the left of the tan 6; = 
Fig. 8 therefore represents hovering instability. As 


—1/n curve of 


§; is increased, the motion becomes gradually very 
stable. This is illustrated in the stability diagrams, 
Figs. 9a, 7, 9b, drawn for tan 6; = — 1/3/12, O, and 
It is seen that the stability for 
6; = O is already adequate, and little can be gained by 


\/3/4, respectively. 
use of positive cocking angles. The figures also indi- 
cate that for nm = +/3 the hovering natural frequencies 
in the three cases are Q) = 0, '/2w, and w, respectively. 
The general variation of /w with 6; at 1 = Ois shown 
in Fig. 8.* 

The stiffening effect of a positive 63 raises the hovering 
natural frequency of the blade and therefore diminishes 
the available percentage critical damping, 
per cent critical = 100nw/2Q = 

50n/V1 + n tan 6; (53) 
This is shown in Fig. 10. The number of vibrations that 
take place before a transient reduces to a certain frac- 
tion of its initial amplitude is therefore increased when 





* In Fig. 8 the curves should be labeled 9/w. 


Errata 


The following errors appeared in a Letter to the Editor by Ralph H. Upson (Journal of the Aeronautical Sciences, 


Vol. 14, No. 8, August, 1947). 
Page 470, first column, first line of the first footnote: 
| “in Mr. Szczeniowski’s Letter.”’ 


Page 470, first column, fifth line from the bottom of the main text: instead of ‘“E, = (v?/2)(mo — mp),’’ it should 
read ‘‘E; = (v?/2)(mo — m).” 
Page 470, top of second column: equation should read ‘““V = at = v/[1 + (g/a)] log. (mo/m).”’ 
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instead of ‘‘in Mr. Krzywoblocki’s Letter,’’ it should read | 
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63 is increased, but the percentage reduction of the tran- 
stent per revolution, 100[1 — exp. (—Mapp.)], is not 
affected by 63. This again shows the relative irrele- 
vancy of the (positive) 6; setting on the transient re- 
sponse of the rotor blade. 

A better appraisal of the effect of the 6; hinge can be 
obtained when one considers the steady-state response 
of the blade to the forced excitation, E(wf). It is to 
be noted that the blade pitch angle #3’ is smaller for 
positive 6; than the manually set value 3, and therefore 
the collective pitch setting must be increased to restore 
the rotor thrust to the value that will support the weight 
of the helicopter. If this is done, then the amplitude 
of the forced response, 8, exceeds (for positive 63) the 
response for 6; = 0. This was shown by the writer's 
colleague, A. Gail, in a 1944 McDonnell Aircraft Cor- 
poration report. Gail also included cyclic pitch varia- 
tion terms in his study. Thus, on the basis of flapping 
response (negative 63's destabilize the transient; large 
positive 63;’s cause undue amplification of the steady- 
state response), a 63; angle of zero or slight positive 
value provides the most desirable hinge setting. 
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Compressive Loading’ 
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SUMMARY 


Use of heavy tapered skins is desirable to meet efficiently the 
load-carrying requirements of high-speed thin-wing aircraft. 
In addition, these skins must be nonbuckling from aerodynamic 
considerations. Consequently, the design of a nonbuckling 
tapered plate acting under various compressive loadings is con- 
sidered. 

Conditions for an efficient plate taper are assumed from physi- 
cal considerations, and the compression instability coefficient for 
a rectangular tapered plate of any aspect ratio is then obtained. 
Design of the plate is then considered for both the elastic case 
and one in which a region is operating at stresses above the pro- 


portional limit. 


SYMBOLS 


length of plate, in 


b = width of plate, in. 

D = flexural rigidity of plate, lbs. per in. 

E = modulus of elasticity, lbs. per sq.in. 

E, = secant modulus, lbs. per sq.in. 

E = areduced modulus, lbs. per sq.in. 

k = instability coefficient in compression 

m = number of waves in a buckled plate 

n = shape parameter of stress-strain curve (reference 2) 

Nz = load per chordwise inch, lbs. per in. 

N, = transverse internal force per inch required for plate 
equilibrium, lbs. per in. 

N,, = shear forces required for equilibrium, Ibs. per in. 

ag = thickness of an elastic plate, in 

t = thickness of an inelastic plate, in. 

w = plate deflection out of its own plane, in 

x = coordinate in spanwise direction, in. 

y = coordinate in chordwise direction, in. 

o* = compressive stress at instability for an elastic ma- 
terial, lbs. per sq.in. 

a = compressive stress at instability for an inelastic ma- 
terial, lbs. per sq.in. 

oo.7 = secant yield stress, lbs. per sq.in. (reference 2) 

8 = loading variation parameter, 8 = In * 

sve 

” = inelastic compression instability parameter, 7 = 
E,/E 

v = Poisson’s ratio 

NoTE: Subscript x refers to the symbol at station x—i.e., 

Diss bay Cas 
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INTRODUCTION 


r THE DESIGN of high-speed, high-performance air- 
craft that employ thin wings, large loads are carried 
at the root section by a wing surface that is expected 
to remain aerodynamically smooth. To attain struc- 
tural efficiency within the limited geometry of a thin 
wing, it is desirable to have the skin act as the main 
load-carrying member. Consequently, it may be nec- 
essary to use wing skins on the order of '/2-in. thickness 
to carry the loads without exhibiting structural in- 
stability. 

In general, instability of the thick skins is tanta- 
mount to failure of the wing, since the skin is the main 
load-carrying member. The large loads on the sup- 
porting structure which are associated with instability 
of the skin are sufficient to cause complete collapse of 
the structure at the onset of instability. Thus, the 
compression allowable of the wing structure becomes 
essentially the critical compressive loading of the weak- 
est section. 

As the loading on the skin decreases outboard of the 
root section, efficient use of material requires that the 
skin be tapered in thickness to just carry the loading. 
In addition, the compressive stress of the skin is a func- 
tion of the thickness of the skin and, consequently, 
varies as the plate tapers toward the wing tip. Thus, 
efficient design of the wing surface resolves into deter- 
mining the correct taper of the skin to meet the loading 
requirement and achieve instability over the entire sur- 


face. 
EFFICIENT TAPER 


As a consequence of the preceding discussion, it is 
desired to design, in an efficient manner, a tapered 
plate that is acting under the loading distribution 
shown in Fig. 1. By efficient is meant a plate that will 
tend to possess minimum weight under the loading dis 
tribution of Fig. 1 and at the same time retain struc 
tural stability under that loading. 

In general, when a tapered plate has attained the 
state of unstable equilibrium, instability is character 
ized by sensible deflections out of the plane of the plate 
The other portions of the plate te 


in one region only. 
This condi 


main essentially free of such deflections. 
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Fic. 1. Loads acting upon tapered plate. 


tion of instability constitutes an inefficient design, since 
the same loading distribution presumably could be sus- 
tained by a lighter plate tapered in such a manner that 
instability under the specified loading will be charac- 
terized by sensible deflections throughout the entire 
plate. 

Therefore, an efficient taper tending toward mini- 
mum weight for the plate shown in Fig. 1 is defined as 
that distribution of plate material whose instability 
under the specified loading distribution is characterized 
by the achievement of simultaneous sensible deflections 
out of the plane of the plate throughout the entire 
plate. 

The difficulties inherent in a mathematical solution 
for the taper of minimum weight are beyond the scope 
of existing instability analysis. Consequently, assump- 
tions must be made as to a reasonable criterion for an ef- 
ficient tapered plate, and the following discussion is 
offered in support of the assumptions involved in select- 
ing this criterion. 

Instability of a flat plate under uniform compressive 
loading is characterized by a uniform sinusoidal buckle 
pattern. In this case, the critical loading of the plate is 
given by 

Nor, =: (9? D/b*)k’ (1) 
where the instability coefficient, k’, is a constant and 
some function of the plate geometry and boundary 


The flexural rigidity, D, is defined in the 
usual manner as 


conditions. 


D = Et*/12(1 — v?) (2) 


Since an efficiently tapered plate was taken to be 
one that has an everywhere-sensible instability pattern 
under the specified loading, the assumption is now 
made that the distribution of plate material to meet 
this condition is such that the loading at each section at 
the onset of instability is given by an equation of the 
type of Eq. (1). Thus 


N, = (x?D,/b?)k (3) 
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where & is a constant for the plate and some function 
of the plate geometry, the boundary conditions, and, in 
addition, the loading distribution. The rigidity at 
each section is similarly defined as 


D, = Et,?/12(1 — vy?) (2a) 


Since the only variables in Eq. (3) are the loading 
and rigidity, and since an everywhere-sensible buckle 
pattern was attained in Eq. (1), where JN, is directly 
proportional to D, the equivalent condition for an ef- 
ficiently tapered plate probably will exist if this rela- 
tionship between loading and rigidity is maintained. 
Consequently, following the assumption involved in 
writing Eq. (3), the criterion for efficient design is the 
distribution of plate material in which the flexural 
rigidity is section-wise proportional to the specified 


loading. Thus, 


N; No _ D, Do (4) 


LOADING VARIATION PARAMETER 


An aircraft wing under air loads possesses a monotonic 
decreasing bending moment along the span. Conse- 
quently, the compression surface must likewise possess 
a monotonic decreasing loading characteristic, since it 
is assumed to be the main member resisting bending in 
the region of the root. 

An exponential function satisfies those general re- 
quirements and was chosen because it may readily be 


adopted to most known loading distributions. The 
loading is represented by 

, ii te — 2—Bbz/ 

N,/No = £(x/a) = e-*7*\ (5) 


In No/N, f 


B 
Similarly, the rigidity of the plate by virtue of Eq. 
(4) may be written 


D, = Dye~®*/4 (6) 


and the required plate thickness is given by 


fs / / -_ 
fo _ toe (8 3) a) (7) 


Determination of the instability coefficient for a 
tapered plate whose flexural rigidity is proportional to 
the required loading must now be considered. 


THE COMPRESSION INSTABILITY COEFFICIENT 


The forces in the tapered plate shown in Fig. | must 
satisfy the following equilibrium equations: 


(ON,/dx) + (ON;,,/d,) = 0) 


, Py (3) 
(ON,,/Ox) + (ON,/0,) 


Il 
my 


The solution of Eqs. (8) is sought for the following 
boundary conditions: 


Ni(y =0)=0 ) 


a ‘ (9) 
N,(y = +5/2) = Of 
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Fic. 2. Compression instability coefficient of an exponentially tapered plate whose rigidity is proportional to the loading. 
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Fic. 3. Elastic thickness required for a given compressive loading. Eq. 
: ( 
These boundary conditions are for a single-cell canti- 1 al teats . (ow? . z Ow Ow No: 

reese Anne ; = N.({—]} + 2Ney (| — + 
evered box beam in bending. 29) 0 J -we2 ox Ox OY 3 
The resulting shears and transverse forces required ow\? 1 {2% (2 {/dw . dw? (2 
for equilibrium, using Eq. (7), are N, |- ) Jaw dy = | / D — +- *) - = 

g Eq. (7 ”\dy 2S 0 S 12 7W\dx? © dy’ 
- + / _B( i 
Ni = No(8/a)ye B(z/a) O°w O2w ow \2 \ 
2/b2 oy (10) i") tease ~ ae e 
Ny a Me (f) E = ) -B(2/a) Ox? Oy? Ox Oy ) 
a 3 Z zs , ; ee 
Following the well-known Rayleigh Principle, the | 
The energy equation to be solved for the determina- use of an assumed deflection form that satisfies the re- | ; 

tion of the compression instability coefficient is given quired boundary conditions for the deflection of the | Simy 
plate will yield a reasonable estimate of the instability coeff 
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DESIGN OF A TAPERED 
coefficient. Since the plate is taken to be simply sup- 
ported along all four sides, the following form is assumed 
for w: 


(12) 


W = Wo sin (mmrx/a) cos (ry/b) 


Eq. (12) is the exact deflection form for a simply 
supported, rectangular plate of constant thickness 
under uniform compressive loading acting in the plane 
of the plate. In the event that large shear forces exist 
in the plane of the plate, comparable in magnitude to the 
applied compressive forces, the assumed deflection is 
no longer an accurate description of the plate deflection. 
It will be noted, however, that in the practical case the 
loading variation parameter [8 = In (0/N,)] is usually 
small and that the plate aspect ratio (a/b) is usually 
large, in which case examination of Eq. (10) will show 
that the shears and transverse forces tend to approach 
zero and Eq. (12) is again a reasonable deflection form 
for the plate at instability. 

Employing Eq. (12), the following integrals are re- 
quired for solution of the energy equation, Eq. (11): 


/20 sin? (axy/b) 
pe (ry 2 te be 
_b/2 cos~* (ry/b): ‘ 
“a | 
SR y sin (2ry/b)dy = b?/2r 
b/2 (72 2 5. Sat 
7. ee b3 [x 
aa af sin? —dy = —{|— — 1 | 
Si e-P*/"dx = af(1 — e~*)/6] | 


. 2mm 
I e-P2/4 ogg dx = 
0 a 
, 1—e~ | (8/2rm)?. | 
B 1 + (8/2xm)? 


a 9 - 
, . 2mrx 
I e-Pz/a cin = dx = 
0 a 
B/2rm 


1 —e* | 
ae — ———_____—_ 
B E + (8 ney J 


Substituting Eqs. (5), (6), (10), (12), and (13) in 
Eq. (11), the energy equation becomes: 


ab 1 — eT /mr\? 
ae: = 

< B ( ) + 
(2) “ (2): ie 1)(£) 

2a 2\a 4\3 a! lop ab1 —e~* 


1+ (6/2xm)? 


(ey + (Tea oO) 


1 + (8/2rm)? 














| 


= 





(14) 





Simplifying Eq. (14) and solving for the instability 
coefficient, k, where k = Nob?/ Dor? 
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a \* mb\? |? B\? 
a) | + 0-45) 
mb a mn 


Ee 
m/ \12 4n° 


The variation of k as a function of a/b for various 8 
values is presented on Fig. 2. This instability coef- 
ficient is used in the expression for the compressive 
loading 


N, = w°Ekt,3/12(1 — v?)b? (3a) 


and is an overall coefficient for the entire plate. It 
should be noted that for 8 = 0 the solution corresponds 
to that for a flat plate. 


DETERMINATION OF TAPER 


Elastic Plate 
From Eq. (7), the required section thickness to ef- 
ficiently sustain a loading, N,, at station x is given by 
t,* ai to* e—(8/3) @/a) (7) 
and utilizing Eq. (3a) 
to* = (Nob?12(1 — v®)/w*kE)”3 (16) 


Eq. (16) is plotted in Fig. 2. 
Similarly, the stress at each section of the plate at 
the onset of instability is given by 


o;* = N,/t,* (17) 
or 
o;* = ape 28/8) (2/0) (18) 
where 
oo* = No/to* (17a) 


Inelastic Plate 


At large loadings for which it is anticipated the plate 
will be designed, it is probable that the operating stress 
at the thick end may exceed the proportional limit of 
the material. In this case, the actual operating stress 
of the plate will be lower than that predicted by elastic 
theory, and it will be necessary to increase the plate 
area to meet the loading requirement. 

The problem of plate stresses exceeding the propor- 
tional limit may be discussed by considering plate 
rigidity. The rigidity of a compressed plate in bending 
is a function not only of the plate thickness but also of 
the compressive stress of the plate. Thus, under 
stresses exceeding the proportional limit of the material, 
the plate will possess a lower rigidity to bending than 
in the elastic region. Consequently, in order to meet 
the loading-rigidity relationship given in Eq. (4), it 
will be necessary to increase the plate thickness in the 
same proportion that the apparent elastic rigidity is 
reduced by exceeding the proportional limit. The re- 
quired condition may be stated as: 
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Fic. 4b. Correction factor for thickness for stresses above the 


proportional limit based on secant modulus correction. 


E(t,*)?/12(1 — v2) = Et,3/12(1 — v2) —(19) 


Here, the usual assumption is made that the entire 
effect of compressive stresses above the proportional 
limit is limited to changes in the effective modulus of 
the material, leaving Poisson’s ratio unchanged. Thus, 
the required thickness may be written in terms of the 
predetermined elastic thickness as: 


t, = t,*(E/E)”3 (20) 


where £ is some reduced modulus and a function of the 
stress level at which the section is operating at the onset 
of instability. 

The determination of the true value of the reduced 
modulus, /, for the combined loadings at which the 
plate section is operating is beyond the scope of existing 
inelastic analysis. Consequently, for the purposes of 
practical solution, the assumption is made that a pos- 
sible value for the reduced modulus for aluminum al- 
loys is given by the secant modulus of the plate ma- 
terial in pure compression.' Thus, 
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t, = t,*(E/E,) = t,*(1/n2)” (20a) 


The true operating stress of the section is given by 
the loading condition that 


N, = ort, = o,*t,* (21) 


or 


or = o;*(nz)® (22) 


Thus, the correction applied to the elastic values of 
both the thickness and the stress required to sustain a 
specified loading is given in terms of the ratio of the 
elastic modulus to the secant modulus. 

It should be noted that, in correcting for the thick- 
ness for stresses beyond the proportional limit by this 
procedure, no change in the value of the instability 
coefficient is introduced by virtue of Eq. (19), which 
left the value of the effective plate rigidity unchanged. 
Thus, with the aid of the material stress-strain curve 
and the solution based on the elastic theory, a correc- 
tion for inelastic effects has been obtained. 
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compressive loading. See Table 1: 75S-T sheet; b = 9./ 1. 
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TABLE 1 


Calculations for Designing an Efficiently Tapered Plate Under an Assumed Compressive Loading 




















———_—_Material 75S-T Plate ——$$<—<—<—<—$ 1 ——————_Tapered Plate Parameters——— ———<—<—<—_ _----—— —Constants———————_ 
E = 10.0 X 106 Ibs. per sq.in. No = 33,400 Ibs. per in. a/b = 3 b2] 12(1 — v2) 
vy = 0.33 Ng = 12,280 Ibs. per in. k = 4.03 (Fig. 2) - = a] 
oo.7 = 70,000 Ibs. per sq.in. B =I1n (No/Nea) = 1 6 = 9.70 in. A = 2.58 XK 1078 
n =10 
1 2 3 4 5 6 7 8 9 10 
Nz = = be Cr, 
x Noe ~B(x/a) or* oz* t,* tz, Lbs. per 
a Lbs. per In. AN; t,* o,* 00.7 (n = 10) (n = 10) in. Sq.In 
Fig. 3 Eq. (17) @) + @) + 
@) xa or (3)'/a Q/@ @)/o0.7 Fig. 4a Fig. 4b @xG ®©®x@ 
0.0 * 33,400 0.08617 0.441 75,740 1.082 0.906 1.104 0.487 68,620 
0.1 30,200 0.07792 0.427 70,720 1.010 0.931 1.074 0.459 65,840 
0.2 27,350 0.07056 0.413 66,220 0.946 0.951 1.051 0.434 62,970 
0.3 24,750 0.06385 0.400 61,870 0.884 0.967 1.034 0.414 59,830 
0.4 22,400 0.05779 0.387 57,880 0.827 0.979 1.021 0.395 56,660 
0.5 20,250 0.05225 0.374 54,140 0.773 0.987 1.013 0.379 53,440 
0.6 18,325 0.04728 0.362 50,620 0.723 0.992 1.008 0.365 50,210 
0.7 16,600 0.04283 0.350 47,430 0.677 0.997 1.003 0.351 47,290 
0.8 15,010 0.03873 0.338 44,410 0.634 0.998 1.002 0.339 44,320 
0.9 13,580 0.03504 0.327 41,520 0.593 0.999 1.001 0.327 41,480 
1.0 12,280 0.03168 0.316 38,860 0.555 1.000 Bd 


38,860 





Analytic Expressions 


The method of Ramberg and Osgood,’ of specifying 
mechanical properties of the plate material by three 
parameters, conveniently provides an analytic expres- 


sion for 
= tin.” 
—-=—=1+4 Pod 
Nr 4s 00-7 


1 3/e,° «.\""* 
—=]+4 2 (% ‘s) 
Nr 4 \G0.7 Gr 
Combining Eqs. (22) and (23), the following analytic 
corrections for critical stress and thickness are obtained. 


“J / 


or 


(23) 


—* o-*(7 o,*\3 1/n—1 
aoa a 
00-7 Or 3 Or , 

o,* te 7 [ te 3 1/n—1 

— = —<{_-|{/—]} —]1 25 
90-7 =e () } ” 


Both these equations are plotted in Fig. 4 as a function 
of o,*/09.7. 


DESIGN PROCEDURE 


The procedure for designing an optimum tapered 
plate for an assumed loading distribution is illustrated 
in Table 1 and Fig. 5. The design for an elastic plate 
is first determined, and then the necessary corrections 
are applied to the region of the plate stressed above the 
proportional limit. 


CONCLUSION 


(1) A design for efficient taper may be obtained by 
a plate whose rigidity at each section is proportional to 
the compressive loading at that section. 


(2) The plate instability coefficient for an efficiently 
tapered plate of any aspect ratio acting under various 
exponential compressive loadings has been determined. 

(3) The theory apparently may be applied to both 
elastic and inelastic materials and used for cases in 
which part of the plate material is operating at stresses 
above the proportional limit. 

(4) For the types of loading usually encountered in 
wing design problems, it has been found from experi- 
ence that a taper with a slight concave curvature is 
more efficient than a linear taper. 
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Discussion 


F. R. Shanley, Division Engineer in charge of Engi- 
neering Research, Lockheed Aircraft Corporation: 
“In extending the analysis to the case of the inelastic 
plate the authors use the secant modulus for the ef- 
fective modulus of elasticity. Although Mr. Gerard 
showed in reference 1 that this method gave results 
that were in good agreement with available test 
data, there is no theoretical basis for using the secant 
modulus in connection with the buckling of plates and 
shells. 

“In a study of this problem made by Lundquist, in 
N.A.C.A. Technical Note No. 686 (‘Local Instability of 
Symmetrical Rectangular Tubes Under Axial Compres- 
sion’), several alternative expressions for the effective 
modulus were considered. The effective modulus will 
have a value that is a function of both the tangent 
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modulus and the elastic modulus, as shown in Timo- 
shenko’s Theory of Elastic Stability, pages 386 and 387. 
For example, one simple expression for the effective 
modulus is 1/ EE, where E, is the tangent modulus. 
(This agrees with Timoshenko’s theory, except that the 
tangent modulus has been substituted for the reduced 


modulus.) Figs. 4a and 4b could have been derived on 


this basis. 
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“The only logical claim that could be made for the se- 
cant modulus is that its value must be somewhere be- 
tween the values for the elastic and tangent moduli, 
This does not seem to be a sound basis for the extension 
of plate buckling theory into the inelastic range. Fig. 5 
of the paper shows that the correction factor can havea 
large effect on the required plate thickness and illus- 
trates the importance of this phase of the problem.” 


Letters to the Editor 


Dear Sir: 

In a Letter to the Editor,! Allen E. Puckett and Ting Yi Li 
point out that Busemann’s approximate solution for isentropic 
Prandtl-Meyer flow 

@ 
Ap =q = C,0" 
n=1 
is incorrect in the third term, and they present the correct ex- 
pression for C;. E. V. Laitone had previously noted that the 
third term, }; = C; — D, of Busemann’s oblique shock wave 
power series was in error and presented the correct value of this 
coefficient in his recent paper.2. Using Laitone’s value of b; and 
their value of C;, Puckett and Li obtained an expression for D 
which, however, appears to be in error. We believe that the 
“correct”’ value should be 
ea worl? s 23( ut faery ad = 
12 (M? — 1)” 4 5 — 37 5 — 3y 
REFERENCES 
1. Puckett, Allen E., and Li, Ting Yi, Letter to the Editor, Journal of the 


Aeronautical Sciences, June, 1947. 
2. Laitone, E. V., Exact and Approximate Solutions in Two-Dimensional 


Oblique Shock Flow, Journal of the Aeronautical Sciences, January, 1947. 
A. KAHANE 
LESTER LEES 
Aeronautical Engineering Laboratory 
Princeton University 


Dear Sir: 
We originally derived the expression for D in the form of 
terms in descending powers of \/, which is actually soniewhat 


neater. This form was correct and is given below. 


(y+ 1)M* [5 — 3y 

CS SS M+ + (y — 3) M?+2 

12(M? — 1)7/ 4 ’ 

In transcribing that equation to a form similar to Busemann’s 

expression we made the error discovered by Kahane and Lees. 

Fortunately, the numerical values given in our Letter to the 

Editor were computed from the original correct form and are 
therefore right. 

ALLEN E. PUCKETT 
Daniel Guggenheim Aeronautical Laboratory 
California Institute of Technology 








Dear Sir: 

On page 540 of the article ‘‘On the Steady Fiow of a Gas 
Through a Tube with Heat Exchange or Chemical Reaction”’ 
(Chambre and Lin, JouRNAL OF THE AERONAUTICAL SCIENCES, 
Vol. 13, No. 10, page 537, October, 1946), it is stated: ‘‘A rapid 
approach to the sonic state can only be achieved if the reaction ab- 
sorbs heat. However, the authors have so far been unable to find 
any practical application of this case.” 


In 1925 the present writer developed the basic principles of 
gaseous reactions with particular reference to engine combustion 
processes, flame-front movements, the rise of detonation, and 
the action of dopes. Many of those principles have now re- 
ceived acceptance. An outline of the theory was published 
(“Heating Fuels for Injection Engines,” by Edward Adams 
Richardson, The Pennsylvania State College Bulletin, The School 
of Engineering, Technical Bull. No. 16, April, 1933). The 
methods therein used are found in part in such modern sources 
as Introduction to Chemical Physics, by Professor Slater of the 
M.I.T., and also Explosions- und Verbrennunsvorgdnge in Gasen, 
by Wilhelm Jost. Suffice it that I am not aware of any experi- 
mental work on combustion or detonation which is not consistent 
with, and explainable by, the principles then laid down. 

In the case of detonation, it was shown that the detonation 
wave is characterized by intermolecular impacts of sufficient 
speed to produce a higher energy (absorbing) level of dissociation, 
permitting much greater energy liberation in the flame front. 
From the flame front, particles are projected in advance thereof 
at speeds approaching the upper limit of molecular velocity. 
The detonation wave velocity is a dynamic equilibrium phe- 
nomena in which the energy released and the energy of dissocia- 
tion at the flame front are in balance. Much energy is liberated 
in the flame front (much more than for ordinary combustion 
waves), and a considerable amount is still liberated behind the 
flame front. 

It will be seen that, so far as the process immediately in ad- 
vance of the detonation wave front is concerned, the chemical re- 
action is characterized by heat absorption. However, the sonic 
velocity is not that of the material ahead of the flame front but 
rather that in, or immediately behind, the flame front; for det- 
onation waves are characterized by supersonic velocities relative 
to the gases ahead of the detonation wave front. 

It is hoped the above remarks may prove helpful in further 
studies of this sort, and they suggest that the change in the 
properties of the fluid between the point immediately in advance 
of the detonation wave and immediately behind may be ex 
tremely large. It is possible the differential equations may not 
permit of reasonable studies in this flame-front region. 

It might be added that the particular fuel process for internal 
combustion engines, based on such principles, was tested by the 
N.A.C.A. and reported upon in “Influence of Fuel-Oil Tempera 
ture on the Combustion in a Prechamber Compression-Ignition 
Engine” (N.A.C.A. T.N. No. 565, April, 1936) and that all of the 
expected benefits were shown to exist. The purpose of the fuel 
process was to permit the use of any fuel, irrespective of its 
chemical character, without detonation and without smoke even 
with air deficiency, while permitting use of the fuel injection 
process without reliance upon the heat of compression of the com-~ 
bustion chamber air for ignition of the fuel. It is regrettable that | 
the importance of the process was not then recognized. 


Epwarp ADAMS RICHARDSON 
Publications 


Bethlehem Steel Company | 
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